There are three fundamental operations that I think any GIS operator should understand (working with **projected coordinates**):

- Calculate the
**distance**between two**points** - Work out the
**heading**(direction, in degrees) from one**point**to another - Calculate a new
**point**at a given**distance**and**heading**from an existing point

In fact, whenever I hired a new member of staff, the first task that I asked them to complete was to create a simple spreadsheet demonstrating that they could apply simple trigonometry to solve these three problems. It was very rare that anyone would be able to do this straight away, but actually it is quite simple - the skills required are nothing more than we learn in high school maths! Most importantly, understanding these few things gives you some great insight into the inner workings of GIS, and you can apply this understanding to solving more complex GIS challenges.

I know that this kind of stuff isn’t everyone’s cup of tea, but with a bit of effort you can get your head round this stuff- and it really does come in handy more often than you would think if you end up doing a job that involves GIS (or any type of mapping whatsoever) - I can’t tell you how many times I have used these!

Here I will explain each of them in turn:

## Distance

You should remember this from **week 2** and **week 3** that when working with projected coordinates, a distance calculation is nothing more than a straightforward application of **Pythagoras’ Theorem**:

where *a* is the length of the *hypotenuse* (the longest side), and *b* and *c* are the other two sides of a right-angled triangle. For for a distance calculation, therefore, this could be better described as:

where \(\Delta\) (delta) means difference (so \(\Delta x\) would be the difference in the x coordinates).

Handily, the Python `math`

library has a convenience function built into it called `hypot`

that calculates the length of the hypotenuse if you give it \(\Delta x\) and \(\Delta y\) as arguments:

```
from math import hypot
def compute_distance(location1, location2):
return hypot(location1[0] - location2[0], location1[1] - location2[1])
```

If you don’t want to use the built in function, you can calculate it yourself using the `sqrt()`

and `pow()`

functions instead:

```
from math import sqrt, pow
def compute_distance(location1, location2):
return sqrt(pow(location1[0] - location2[0],2) + pow(location1[1] - location2[1],2))`
```

Or you can even do it without any imported functions whatsoever:

```
def compute_distance(location1, location2):
return ((location1[0] - location2[0])**2 + (location1[1] - location2[1])**2)**0.5
```

*Of course, if you want to do this for geographical coordinates; I would recommend using the Inverse Vincenty Equation to calculate the distance over an ellipsoid. You can read more about this (as well as various methods for spherical distance calculations) in the distance page; and see a brief example on the handy hints page.*

## Heading

To work out the **heading** (direction) between the two locations, we need to use a little bit of **Trigonometry**. You might remember **SOH CAH TOA** from school, which is a simple pneumonic that helps you remember whether you need to use `sin`

, `cos`

or `tan`

:

**SOH CAH TOA** is intended to help you work out which formula you need to solve a triangle. Each of the three formulas contains a combination of the length of two of the three sides of the triangle, and the angle \(\theta\). All that you have to do is select the one that includes both the information that you have access to and the information that you want.

For example, we are wanting to work out the **direction** between two points (those at either end of the Hypotenuse), which in the above diagram means working out the angle \(\theta\). If we are given the coordinates of the two points, then we know the length of the **Adjacent** side (\(\Delta x\), green in the above image) and the length of the **Opposite** side (\(\Delta y\), red in the above image), because we can simply subtract the coordinates to get the differences between them.

We therefore want to use the **TOA** equation, as this is the only one that includes the **angle** (\(\theta\)), the **Opposite** and the **Adjacent**:

Or, in more geographical terms:

\[Tan(\theta) = \frac{\Delta y}{\Delta x}\]Because we want to know \(\theta\), rather than \(Tan \theta\), we can rearrange this equation like this:

\[\theta = Tan^{-1}(\frac{\Delta y}{\Delta x})\]\(Tan^{-1}\) (**Arc Tangent**) is the opposite of \(Tan\) (**Tangent**), and its functionality is provided in the `atan`

and `atan2`

functions of the Python `math`

library. `atan2`

is just a convenience function that allows you to pass in \(\Delta x\) and \(\Delta y\) as two separate arguments, and it handles the division for you.

Irritatingly for geographers, maths libraries in Python (and almost all other programming languages) calculate their results in **radians** rather than **degrees**, so you need to convert them using the `degrees`

function in the Python `math`

library:

```
degrees(atan2(location1[1] - location2[1], location1[0] - location2[0])
```

Another way in which `atan2`

might confuse a geographer is that it returns the angle measured from **south** on a scale of $-180$ to $180$, where we (as geographers) want the angle from **north** on a scale of $0$ to $360$. Fortunately - this is pretty easy to convert, we simply add $180$ to \(\theta\):

```
180 + degrees(atan2(location1[1] - location2[1], location1[0] - location2[0])
```

Finally, there is also a mathematical trick to play to ensure that the number you end up with is between 0 and 360. To achieve this, you need to add `360`

to your result, and then calculate the **modulus** of your result and `360`

using the **modulo** oparator (`%`

). The **Modulus** is the remainder after you divide the numbers on each side of the operator, so, for example:

```
100 % 40 = 20
```

because `40`

goes into `100`

twice (`80`

), leaving a remainder of `20`

. Here it is in context:

```
(180 + degrees(atan2(location1[1] - location2[1], location1[0] - location2[0])) + 360) % 360
```

To combine all of what we have discussed, then - the magic formula for calculating the heading between two points in **degrees clockwise from north** (which is what geographers expect…) is:

```
from math import degrees, atan2
def compute_heading(location1, location2):
return (180 + degrees(atan2(location1[1] - location2[1], location1[0] - location2[0])) + 360) % 360
```

*If you want to do this for geographical coordinates; I would recommend using the Inverse Vincenty Equation to calculate the forward and back azimuth (as well as distance) over an ellipsoid. You can read more about this in the handy hints page. Remember that the direction is presented as forward azimuth and back azimuth - as the direction relative to north changes throught the trajectory.*

## Point Offset

Another rather common problem in GIS that can be solved by **Trigonometry** is that of the **offset**: calculating a point at a given distance and direction from another point *(e.g. what are the coordinates of a point 3km due east of my house)*.

Fortunately, this is also quite straightforward once you understand **SOH CAH TOA**. Let’s have another quick look at the image:

This time, we have the length of the *Hypotenuse* (red, the distance betwen the two points), and we need to calculate the location of the second point, which is simply:

Our task is then to calculate \(\Delta x\) (**Adjacent**, green) and \(\Delta y\) (**Opposite**, blue) and add them to \(x\) and \(y\) respectively!

Firstly, we need to decide which formula to use. For \(x\), we already have \(\theta\) and the length of the Hypotenuse (**H**) and we need the length of the Adjacent (**A**), so we need the **CAH** equation:

Or, in geographical terms:

\[cos(heading) = \frac{\Delta x}{distance}\]Rearranging this to apply to our problem, we get:

\[\Delta x = cos(heading) \times distance\]So, to get the coordinates of \(x\):

\[x_2 = x_1 + (cos(heading) \times distance)\]As with `atan2`

above, the `cos`

function in the Python `math`

library expects angles to be in **radians**, so we must convert our heading from degrees to radians using the `radians`

function from the Python `math`

library.

```
from math import cos, radians, sin
def compute_offset(location, distance, direction):
return (location[0] + cos(radians(direction)) * distance,
location[1] + sin(radians(direction)) * distance)
```

*If you want to do this for geographical coordinates, I would recommend using the Forward Vincenty Equation to calculate the offset location over an ellipsoid. You can read more about this in the handy hints page.*

## Testing

Let’s see if it works:

Pick two random sets of coordinates (e.g. I tend to use

`345678, 456789`

and`445678, 556789`

) and`computeDistance`

and`computeHeading`

functions.Then calculate the offset point from your first set of coordinates, using the distance and heading values from the above calculation.

If all has gone to plan, the coordinates that you get back should be your second coordinate pair! Here it is in action:

```
from math import hypot, degrees, atan2, sin, cos, radians
def compute_distance(location1, location2):
"""
Compute the distance between two locations using Pythagoras' Theorem
"""
return hypot(location1[0] - location2[0], location1[1] - location2[1])
def compute_heading(location1, location2):
"""
Compute the heading (direction) between two locations using trigonometry
"""
return (180 + degrees(atan2(location1[1] - location2[1], location1[0] - location2[0])) + 360) % 360
def compute_offset(location, distance, direction):
"""
Compute the location of a point at a given distance and direction from a specified location using trigonometry
"""
return (location[0] + cos(radians(direction)) * distance, location[1] + sin(radians(direction)) * distance)
# set locations
location1 = (345678, 456789)
location2 = (445678, 556789)
print(f"\nlocations:\t\t\t{location1}, {location2}")
# calculate distance and heading
dist = compute_distance(location1, location2)
print(f"distance:\t\t\t{dist:.2f}")
head = compute_heading(location1, location2)
print(f"heading:\t\t\t{head}")
# calculate location 2 from location 1 and the distance and heading
new_location = compute_offset(location1, dist, head)
print(f"location 2 (re-calculated):\t{new_location}")
print(f"Expected result:\t\t{location2 == new_location}\n")
```

Which returns:

```
locations: (345678, 456789), (445678, 556789)
distance: 141421.35
heading: 45.0
location 2 (re-calculated): (445678.0, 556789.0)
Expected result: True
```