GIS Algorithms and Formulas for Spatial Analysis

1. Introduction • Previously we looked at computer languages. • Today to give a feel for how programmers think we will look at some formulae and algorithms used in GIS. • An lagorithm can be thought of a set of rules with branches and loops.

2. Measuring Distance • Cartesian distance (Pythagoras): • Geographic coordinates:

3. Intersection Of Straight Lines • Formula of a straight line: y = mx + c • Gradient: m = (y2-y1) / (x2-x1) • c is the intercept • At intersection of two lines: • y = m1.x + c1 = m2.x +c2 • \ (m1-m2).x = (c2-c1) • \ x = (c2-c1) / (m1-m2) • At intersection: • y = m1.(c2-c1) / (m1-m2) + c1 • \y = (m1.c2 - c1.m2) / (m1-m2)

4. Do Polylines Intersect? • We can see whether two lines intersect by looking at them. • Computers require a set of rules to test whether they intersect.

6. Do Polylines Intersect? (2) • Could calculate the intersection of every segment in line 1 with every segment in line 2, but inefficient. • Can speed things up using a minimum bounding box around each polyline. • The boxes are imaginary – i.e. they are not actually drawn.

8. Do Polylines Intersect? (3) • Lines cannot intersect if boxes do not overlap (e.g. lines 1 and 2) • Overlap can be easily tested: If the smaller of the two maximum x values is less than the larger of the two minimum x values, then the two enclosing boxes could not overlap. Likewise for the y co-ordinates. • An overlap does not necessarily mean the lines intersect (e.g. lines 1 and 3) • Must do further test if boxes overlap • Further tests involve testing for overlap between minimum bounding boxes around each segment in the two polylines

10. Do Polylines Intersect? (4) • The box for segment 1 overlaps the box for segment 2, likewise the boxes for segments 2 and 3. However, only the latter intersect. • To test, calculate the intersection point using the formulae above, then test whether this lies inside one of the bounding boxes.

11. Area Of A Polygon • Calculating the area of an irregular shaped polygon is surprisingly simple, although the derivation is by no means obvious. • Can think of the polygon as having a top side and a bottom side. • The area of the polygon is the area between the top side and the x axis minus the area between the bottom side and the x axis.

13. Area Of A Polygon (2) • To calculate the area between a polyline and the x axis, you can think of the area divided into a number of trapeziums:

14. Area Of A Polygon (3) • The area is each trapezium can be calculated by regarding it as a right angled triangle sitting on top of a rectangle: • (x2-x1)*(y2-y1) / 2 plus (x2-x1)*y1 • This simplifies to: • (x2-x1)*(y1+y2) / 2

15. Area Of A Polygon (4) • Deciding whether a segment is part of the top or the bottom of the polygon could be problematic, but the problem resolves itself if you work round the polygon in a clockwise direction. • Doing this, by the time you get to the bottom of the polygon, the first x coordinates will be lower than the second and therefore produce a negative area. • The area of the polygon is therefore simply:

16. Location Of A Centroid • Locating the centroid of a polygon is not very elegant. • Need to construct imaginary triangles between each segment and a point opposite. • Centroid is the weighted average of the centroids of all the triangles.

17. Point In Polygon Test • To determine whether a point lies within a specified polygon, we construct an imaginary line from the point to either the x or y axis. • We then count how many times it intersects the polygon. • If it intersects the polygon an odd number of times, then the point must lie inside; if it does not intersect the polygon at all or intersects it an even number of times the point must lie outside.

19. Point In Polygon Test (2) • To count the number of intersections, we could compare the line with every segment of the polygon using the method discussed earlier. • However, it can be done more simply. All points on the line have the same x coordinate. Any segment on the polygon intersecting this line must have one x coordinate with a lower value and one with a higher value. • We just need to copare each sucessive pair of x cordinates, and count the number of intersections.

20. Polygon Allocation • In many situations we need to decide which of several polygons a point lies in. • We could do the above test with each polygon in turn until we get a hit, but this is not very efficient. • To speed things up we can construct an imaginary mimimum enclosing rectangle around each polygon, and then test whether the point lies within the box. • This simply requires us to compare the cordinates of the point with the minimum and maximum x and y coordinates of the rectangle.

21. Polygon Allocation (2) • Having quickly eliminated the majority of polygons, we could use the point in polygon test above on the remaining candidates. • If the data are topologically encoded, then we can speed this up even further. All we need do is identify the first intersection between our line and a segment of the polgon, and then read the topological information to establish which polygon is to the left (or right) of the arc that segment forms part of.

22. Polygon Overlay

23. Affine Transformations • Affine transformations may be used to convert from one Cartesian system to another (e.g. digitised coordinates to a specific projection). • There are three basic transformations: • Translation • Scaling • Rotation

24. Translation • Translation simply refers to the movement of a point from one location (x,y) to a new location (x',y'). • If the displacement in the x and y directions is Tx and Ty respectively, then the formulae are simply: • x' = x + Tx • y' = y + Ty

25. Scaling • To scale a point (x,y) we need to apply two scale factors Sx and Sx. The formulae are: • x' = x. Sx • y' = y. Sx • Scaling will result in displacement. We therefore need to identify a fixed point (Fx, Fy) which can be used to 'anchor' the feature (e.g. its centroid). • Formulae become: • x' = x. Sx + Fx (1- Sx) • y' = y. Sy + Fy (1- Sy)

26. Rotation • To rotate a point about the origin by an angle q (measured anticlockwise), the formulae for the new point are: • x' = x.cosq - y.sinq • y' = x.sinq + y.cosq • As before, we need to define a fixed point which will act as the pivot point. The full transformation is equivalent to moving this fixed point to the origin, doing the rotation, and then restoring the fixed point back to its original location.

27. Changing A Coordinate System • We need to identify three things: • a common point that can be defined in both coordinate systems, • the scale factors in each dimension, and • the orientation of the axes in one system relative to the other.

28. Conversion Steps • Conversion involves the following steps: • Translate by -Ox, -Oy (where Ox,Oy is the origin of B expressed in A units). This is equivalent to moving the origin of B to the origin of A. • Rotate by -q, where q is the angle between the x axes in the two systems. This brings the axes of B into alignment with those of A. • Scale by Bsx, Bsy (where Bsx is number of x units in A equivalent to 1 x unit in B, and Bsy is the same for the y coordinate units)

29. Vector To Raster Conversion • Procedure is to ‘switch on’ the cells intersected by a line. • However, wish to keep lines as thin as possible without any gaps. • Number of cells required is: • n = Dmax + 1 • where Dmax is the larger of Dx and Dy, defined as: • Dx = abs (x2 - x1) • Dy = abs (y2 - y1)

31. Vector To Raster Conversion (2) • To identify which n cells to 'switch on', we can calculate the size of the steps we need to move along the line in the x and y directions using the formulae: • incx = Dx / Dmax • incy = Dy / Dmax • Polygon fill algorithm could be used to colour a raster polygon.

32. Raster To Vector Conversion • Steps involved: • Thin raster line to minimum necessary – various algorithms (e.g. peeling, distance skeleton) • Identify nodes • Trace line – i.e. replace each cell by a vector vertex • Remove unnecessary vertices