Converting distance/azimuth to a real position

Introduction

This will be a short tutorial on how to take data from a field survey and prepare it for a GIS.

In forestry, field crews record the location of trees using azimuth and distance from a fixed position. The crew sets up a stake at the center of the plot, then they take out a compass and face north. They work through the plot like a clock, labelling each tree with a metal tag. One person will stand at plot center with the compass and measure the angle each tree is from North, referred to as the azimuth. Another person will take a logger's tape and record the distance to each tree from the stake. This is the standard way to do fieldwork in the forest, because it's inexpensive, easy to replicate, and the margin of error is adequate for future fieldwork.

Concept

I will show a little tutorial on how to take this data acquired in the field and prepare it for analysis. I will use R and the sf library. Firstly I will cover the concepts, then show it in code.

We have the following information:

• Absolute position of plot center (i.e. latitude/longitude)
• Relative distance of tree(s) to plot center
• Relative azimuth of tree(s) to plot center

Our task is to get the absolute positions of each tree in the same coordinate system as the plot center.

A typical entry looks something like this

We have an orientation (from azimuth) and a magnitude (from distance), so we have a vector. In order to get the 2D coordinates, we use trigonometry to get the two components of this vector. However, we can't just use our azimuth by itself. The azimuth is measuring clockwise, starting from north. The angles in the unit circle progress anticlockwise, starting from east. Therefore we need to multiply our azimuth by -1 (to go from clockwise to anticlockwise) and then add 90 degrees (to rotate from east to north).

I occasionally see people make this mistake and it's quite frustrating to debug, so hopefully I have helped some of you out! For clarity I will always refer to azimuth (α) as the clockwise angle from north, and theta (θ) as the anticlockwise angle from east.

Common error #1
Forgetting to convert the azimuth α into θ before using a trigonometry function.
Solution: Treat α and θ separately. Check the output of trigonometry functions. When θ=0, cos(θ)=1 and sin(θ)=0.

We now have the positions of all the trees relative to the plot center, using an origin of (0,0). In order to calculate the real world position, we just need to translate the relative coordinates to use the real world origin. In order to do this, we add the x and y coordinates of the plot center to our equation. Since the plot center is the same for all trees, we can do it at the end.

An important caveat is that all these measurements need to use the same units and axes. The default settings of most GPS receivers is in decimal degrees, whereas your measurements will be in meters. Any GIS will have functions to do this transformation for you, but it's good to be aware. I have also seen some people have trouble from reversing the x and y coordinates (latitude is the y axis, longitude is the x axis).

Common error #2
Improperly mixing units of measure.
Solution: Check the CRS of all data before you start, and transform them into the same CRS for processing.

Common error #3
Reversing axes
Solution: Read the headers of all data before you start. Use column names to access your data rather than indexing.

Code

Now that I explained the concept, I will show some R code using the sf library to achieve this. sf stands for simple features and it's a very nice library for working with geospatial data.

I have an example dataframe with some values:

sf provides many different objects to represent geospatial data. In this tutorial we will be using the simplest types: point and multipoint.

First, let's create the plot center as a point. To keep this readable and explicit I won't nest functions. Let's say our plot center is (100, 250).

xy = c(100, 250)
point_xy = st_point(xy)
plot_center = st_geometry(point_xy)

We have now created the geometry we need for our plot center. However, there is one important thing we still need to do: set the CRS (coordinate reference system). The two typical ways to represent a CRS are EPSG codes or PROJ.4 strings. Not every CRS will have an EPSG code, but so many do that EPSG code has overtaken PROJ.4 as the preferred way to describe a CRS. The EPSG code is a convenient representation of the corresponding PROJ.4 string. You can look up EPSG codes on a website such as EPSG.IO. It also includes the descriptions in formats such as PROJ.4, WKT, etc. This can be useful when working with older data, which is quite common in GIS.

To set the CRS code, we use the st_crs() function. Using an EPSG code:

st_crs(plot_center) = st_crs(3310)

Using the PROJ.4 string:

st_crs(plot_center) = st_crs("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ")

I know this may sound like overkill to go into this much detail, but coordinate systems are a very common issue when working with geospatial programming. Now that we have this location, we're ready to translate our trees. We can use st_coordinates() to access the position of our plot center and store them as constants.

plot_x = st_coordinates(plot_center)[1]
plot_y = st_coordinates(plot_center)[2]

Now we apply the function I mentioned above. First, we need to transform our azimuth measures, then convert them to radians for R's trigonometry functions.

trees$theta = trees$AZ * -1 + 90
tree_x = plot_x + trees$DIST * cos(trees$theta*pi/180)
tree_y = plot_y + trees$DIST * sin(trees$theta*pi/180)

To create the multipoint object, we need to provide a 2xn matrix with the coordinate, and finish up by adding the CRS.

tree_mat = as.matrix(cbind(tree_x, tree_y))
tree_points = st_multipoint(tree_mat)
trees = st_geometry(tree_points)
st_crs(trees) = st_crs(3310)

Finally, we export with st_write(). Most GIS data is in the form of shapefiles, but I highly recommend you use a more modern format, personally I use GeoJSON.

st_write(trees, "/path/to/your/outputs.geojson")

And now we're all done! We can do a final sanity check by plotting the trees (green triangle), the plot center (white circle), and the plot boundary (big pink circle).

Conclusion

Hopefully you found this interesting and useful. This is applicable to any data that relies on relative distance/bearing to a known position.