DEV Community: theohartsook

Converting distance/azimuth to a real position

theohartsook — Tue, 28 Dec 2021 03:03:47 +0000

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.

Things I learned to appreciate in graduate school: abstraction

theohartsook — Thu, 02 Dec 2021 01:23:43 +0000

I finished my masters degree over the summer and started a PhD program right after. I have always felt a little bit of the odd one out because my background is in environmental science and GIS rather than computer science or math. However, programming was a big part of my masters and will continue to be in my PhD. I also dabbled with it in my career before going back to school. I would like to share some of what I learned over that time, because I feel that some concepts I struggled to learn in class really began to make sense once I started writing semi-serious code.

I never really got the point of abstraction before. I saw the obvious benefit of being able to use print() and other handy functions, but I never really thought about while writing scripts. Once I started working on a LiDAR processing pipeline, I really struggled. The pipelines are long and complicated, with some steps requiring GUIs, some steps easy to automate, and many steps that need to be done manually first to find good parameters, which can then be automated.

This took me one semester to implement, and another semester to get comfortable with. It took a lot of meetings with advisors, mentors, and kind strangers, but eventually I got us a working pipeline. Then I needed to do the tough job of explaining how and why to use this pipeline.

In my experience LiDAR is really tricky to get started with for people of all kinds of backgrounds. For people in the natural science domain, such as me, it's much more math and computer science heavy than we are used to. For people in the computer science domains, it seems to be the geodetic elements. On top of all this, there is the logistical challenge of file formats, vendor specific software, identifying noise, and so on. There is a lot to explain and it is easy to get confused.

In order to explain a pipeline like this, I really had to start using abstraction to describe it at different levels. This made me think a lot about the difference between the general overview and the implementation level.

This is an example of the high level overview I would give.

And for this example, I zoomed in on specific step (getting the data from the scanner) and included some of the decisions that need to be made.

This requires a lot of decision making and experience really makes a difference here. The process can be automated, but it requires a lot of parameters to be decided. These parameters vary depending on the type and platform of scanner, wavelength of the laser, the weather, the number of scanning positions... It takes a lot of work to get centimeter or millimeter level accuracy!

Obviously this is a lot to explain to people, and it's a lot of details that people don't necessarily need to know right away. Especially if someone is just curious about what LiDAR can do and if it would be applicable for their work. This pipeline can be described at multiple levels of abstraction, so it's up to me to pick the best one for my audience at that time. I often want to go into the details, but "I used PDAL to prepare these rasters" is appropriate.

Debugging in CS50

theohartsook — Fri, 14 Aug 2020 16:25:28 +0000

Over the summer I have been working on Harvard's CS50 EdX course and it's been a really fun experience (working on C++ is on pause until I finish CS50). After working with pointers in C I finally feel like I understand (big picture) how memory works.
One of the assignments is to create a simple box blur. One of the neat features of CS50 is that they provide tests for the code to see if we completed the assignment correctly or not. I really appreciate this as I can't visually assess if a blurred pixel is the right kind of blurred.
I had been pretty dispirited this week because I wasn't sure why my blur was failing. I was getting most of the blur correct, but it would fail for one pixel in one test case and two pixels in the second test case.

I knew I had solved it in a suboptimal way because I was treating each corner and edge differently, rather than working with a buffer or some of the other solutions I saw online later. It's a pretty janky solution, but it's one I came up with independently. I figured there would be some kind of problem here, and my intuition was ultimately right.
My first thought was that since it only seemed to the last row that failed, I should check that I wasn't accidentally taking the average of the wrong color channel or indexing wrong. Those turned out to be correct, no matter how many times I reread them.
And then this morning, when I had a free 15 minutes, I decided I would take a look at it again. I copied the expected output and my output into a table so I could look for exactly what was wrong. This instantly revealed the problem. My corner pixels were correct, but the bottom edge pixels were failing. I checked that section of my code again and everything looked correct... until I finally noticed the mistake.
// bottom edge else if (top_edge == 1) { ... // calculate the average for each color channel }
I was checking the wrong condition. I corrected that, recompiled, reran the tests, and passed. I'm glad I finally fixed that mistake.