Last deck told you which features are related to which. This one puts that to work. A spatial join takes the attributes of one layer and attaches them to another based on where things are, rather than on a shared identifier. That is the operation that turns spatial data into something you can regress. You have county-level nitrogen use and watershed-level water quality with no common key between them, only geography, and a spatial join is what marries the two. We work through every combination of points and polygons, then two refinements: doing the join and the summary in one step, and weighting by how much area actually overlaps. The learning objective on screen is therefore broader than memorising one function: you should be able to overlay two s-f objects, identify the relevant source features for each target feature, and bring their attributes together. The table of contents gives you the route. We start with the spatial-join idea, work through points against polygons, polygons against points, and polygons against polygons, then change the topological relation, combine a join and summary with aggregate, and finish with a cropping join. If the plotting or data-manipulation syntax is unfamiliar, use the ggplot-two and d-plyr primer links on the right before continuing.
The objective of this chapter is to learn spatial operations that involve two sf objects. Specifically,
sf object on another sf object to extract (or join) values from the sf objectSame arrangement as always. One habit that pays off particularly on this deck: after every join, print the result and count the rows. Joins change the shape of your data in ways that are easy to miss, and a join that silently multiplied your dataset fivefold will quietly ruin everything downstream. Several slides here exist precisely to make you look at row counts. Take them seriously rather than skimming to the next piece of syntax, because knowing what a join did to your data is most of the skill. The other thing worth doing is reading the row counts against the counts you started with, every time. The controls on screen let you review actively. Open the three-line menu at the lower left to jump through the table of contents, or press the letter o to see the whole deck in overview. The pale-blue area is an editable code area. Run Code evaluates the entire area, while highlighting only the lines you want and pressing Command-Enter on a Mac or Control-Enter on Windows evaluates just that selection. The two-sheets icon copies the code so you can paste it into R on your own computer. If your edits leave a cell in a confusing state, the reload icon immediately to the left of copy restores the original code. Use those controls to test the row-count habit rather than only reading about it. One more, at the right end of the toolbar: the eye icon hides that cell’s output and a second click brings it back. Some of these results run long, and once you have read one it is just pushing the rest of the slide out of view.
Click on the three horizontally stacked lines at the bottom left corner of the slide, then you will see table of contents, and you can jump to the section you want
Hit letter “o” on your keyboard and you will have a panel view of all the slides
Read the two definitions carefully, because target and source are the vocabulary for the whole deck and they are easy to mix up. The target is the layer you want to end up with, the one you are adding columns to. The source is where the values come from. The little coloured words in the callouts are the hint: values go to the target, from the source. And note the third bullet in the list, about summarising if necessary. That is the part people forget. When several source features match one target feature, you have to decide what to do with them, and that decision is yours.
Spatial join involves all of the following:
Definitions: target layer
The sf layer that has sfgs to which you would like to assign the values of attributes from the source layer.
Definitions: source layer
The sf layer that has sfgs from which you would like to assign the values of its attributes to the target layer.
Four combinations exist, since either layer can be vector or raster, but you only need two of them. Vector against vector is this deck. Vector against raster, meaning extracting raster values to points or polygons, is the next deck and is enormously common in practice, since weather and soil data are rasters. The other two have a raster as the target, and the honest reason we skip them is that your unit of analysis is essentially never a grid cell. It is a field, a county, a farm, a well. Those are vectors.
We can classify spatial join into four categories by the type of the underlying spatial objects:
Among the four, our focus here is the first case (the second case will be discussed later).
We will not cover the third and fourth cases in this course because it is almost always the case that our target data is a vector data (e.g., city or farm fields as points, political boundaries as polygons, etc).
Within vector against vector there are nine combinations of points, lines and polygons, and we cover three. Lines get dropped, and the reasoning is worth hearing: lines are rarely the thing you observe and rarely the thing you extract values from. A road is context, not an observation. That leaves points into polygons, polygons into points, and polygons into polygons, and the next three sections take one each. All three use the same function. What changes is what comes back — one row per target feature in the first case, and potentially many in the other two — and that difference is really the substance of this deck.
As noted earlier, we will look at vector-vector interactions in this lecture.
This category can be further broken down into different sub categories depending on the type of spatial object (point, line, and polygon).
Here, we will ignore any spatial joins that involve lines. This is because objects represented by lines are rarely observation units in our analysis nor the source data from which we will extract values.
Here is the list of the types of spatial joins we will learn.
The first and simplest case. In this example, every well matches exactly one county, so every well gets exactly one set of values, and the join adds columns without adding rows. Point-to-polygon joins do not guarantee one row per point: an unmatched point is retained with missing source attributes by default, while a point intersecting overlapping polygons is duplicated. The function is st_join, target first and source second, and that argument order is the whole interface. As with subsetting last deck, the default relation is intersects, and the polygons that intersect each point determine its matches.
What?
For each of the observations (points) in the target points data,
How?
In order to achieve this, we can use the st_join() function, whose syntax is as follows:
Note
Similar to spatial sub-setting, the default topological relation is st_intersects()
Nebraska irrigation wells as the target, county boundaries as the source. The map shows the setup: thousands of points scattered over ninety-odd polygons. The question this join answers is which county each well is in, which sounds trivial and is exactly the sort of thing you would otherwise do by hand or get wrong. Notice you are not told which county any well belongs to anywhere in the data. That information exists only in the coordinates, which is exactly the point of a spatial join. Any other kind of join would need a shared column, and here there simply is not one to use. Read the plotting code from top to bottom. An empty ggplot call starts the canvas, the first geom-s-f draws ne-counties, and the second draws wells-ne on top so the points are not hidden by the polygons. Theme-void removes axes, ticks, and grid lines because those decorations do not help you see containment. The cell runs automatically when the slide opens, so the displayed map is the visual statement of the join problem. Move to Code to turn that spatial relationship into columns.
We use wells_ne (points) and ne_counties (polygons) data for illustration. Here is the map:
One line. Target first, source second, and the wells come back carrying their county’s attributes. Worth saying explicitly what did not happen: the number of rows did not change. You had one row per well before and you have one row per well after, just with more columns. Hold on to that, because the very next section does something different and the contrast is the thing to learn. If this were all a spatial join ever did, it would be a much less interesting and much less dangerous operation. Print the object and confirm the row count for yourself before moving on.
For each of the points (wells) in wells_ne, the code below will find the polygon (county) in which the point (well) is located, and attach the value of the variables of the polygon (county) to the point (well).
Print the result and you will see statefp, countyfp and name appended to each well. Then do the check on this slide, which is the habit worth forming: filter to one county according to the join, pull that county’s boundary separately, and plot the two together. If the join worked, every point lands inside the outline. If it did not, you will see it instantly. Verifying a spatial join visually costs about four lines and catches mistakes that are otherwise invisible. A mismatched coordinate system stops the join with an error, while correctly aligned layers with no spatial matches return missing source attributes in the default left join. The first filter keeps joined rows whose appended name equals Adams, so those are wells classified as being in Adams County. The second filter applies the same name test to ne-counties to obtain the reference boundary independently. In the verification plot, the county geometry is drawn first and the selected wells second, then theme-void removes non-spatial decoration. That order makes any point outside the Adams outline easy to spot. This check tests the spatial assignment itself, not merely whether the expected columns exist.
Evaluate wells_joined_with_county and you will see that statefp, countyfp, and name variables are appended.
Let’s check if the two datasets are indeed joined based on their spatial locations.
Visualize:
Now turn it around. Counties are the target, wells the source, and the difference is not cosmetic. A point sits in one polygon, but a polygon can contain any number of points, so this direction is one-to-many. Same function, same argument order rule, target first. What changes is what you get back, and the Inspect tab is where that becomes clear. Note the goal stated on the Data tab: average groundwater extraction by county. Getting there needs a second step after the join, which is the last tab of this section.
What?
For each of the observations (polygons) in the target data,
How?
In order to achieve this, we can use the st_join() function, whose syntax is as follows:
Note
Similar to spatial sub-setting, the default topological relation is st_intersects()
The same two layers, swapped round. The extraction figures are invented, generated with runif, since the real ones are not public, but they behave like the real thing for our purposes. What matters is the goal: an average per county. Right now that number does not exist anywhere. Extraction is recorded per well, counties have no extraction figure at all, and the only thing linking them is position. Two steps get you there. The join is the first, attaching every well’s extraction figure to the county it sits in, and the summarise is the second. The join is spatial, while the summarisation is an ordinary non-spatial calculation. The first automatically run cell repeats the map, with county polygons drawn before well points and theme-void clearing the axes. The next cell changes the data. Mutate adds a column called g-w-extraction to wells-ne; n-row supplies one draw for every well, and run-if draws each fake value between the named minimum of fifty and maximum of four hundred acre-feet. The dot refers to the piped wells-ne object when its row count is calculated. The assignment saves the new column, while the outer parentheses also print the updated object so you can inspect it. Move to Code with the target, source, and outcome now clearly defined.
We use ne_counties (polygons) and wells_ne (points) data for illustration. Here is the map:
We create a fake variable that represents groundwater extraction (acre-feet) from the aquifer.
Our goal is to find average groundwater extraction by county.
Identical syntax to before, arguments swapped. That is the whole change. It is worth appreciating how little the interface asks of you: the same function handles both directions and works out what to do from the geometries themselves. But do not let the sameness of the code fool you into expecting the same kind of result. Print what comes back before reading on, and look at how many rows it has compared with the ninety-three counties you started with. If that number surprises you, the next tab explains exactly why, and it is the thing most worth understanding in this deck.
For each of the polygons (counties) in ne_counties, the code below will find all the points (wells) that are located inside the county, and attach the value of the variables of the points (wells) to the polygon (county).
Here is the one-to-many consequence, and it is the most important slide in the section. Each county now appears once per well inside it. Loup county, the first row, has no wells, so it appears once with missing values. Adams county has fifty-three wells, so it appears fifty-three times, with the identical county geometry repeated in every one of those rows. The identical check at the bottom proves that. This is normal and expected for a one-to-many join, but if you did not know to look, you would think your data had been corrupted. The printed joined object has 1,052 rows rather than the 93 county rows you began with, which is the multiplication you are diagnosing. Filter then uses the county FIPS code string zero-zero-one to isolate Adams and saves those 53 rows as adams-county; the outer parentheses print them as part of the same action. Finally, selecting row one and row two with only the geometry column gives two one-row s-f objects, and identical returns true because both carry the same Adams multipolygon. Once you understand why that repetition occurs, move to Post-join processing to collapse the repeated rows to a county summary.
Evaluate county_joined_with_wells and you will see that wellid and gw_extraction variables are appended.
One thing that is different from the previous case is that
For each of the polygons (counties), the resulting dataset has as many observations as the number of wells that intersect with the polygon (county).
If a polygon has no wells inside, then you will simply have a single row of data for that polygon.
For example, the first row is countyfp 115 (Loup county), which has no wells inside it. So we get a single row with wellid and gw_extraction missing. But countyfp 001 (Adams county) has 53 wells inside it, and so appears 53 times.
All the rows there have exactly the same geometry, which is the MULTIPOLYGON that represents the boundary of Adams county.
The second step: group by county and summarise, which is ordinary dplyr with nothing spatial about it. Now read the two bullets, because they explain something you will hit constantly. Summarising an sf object is slow, and the reason is that geometry travels with the data, so summarise unions the geometries of each group as well as averaging the numbers. That is real work you did not ask for. Uncomment the line that drops the geometry and it becomes a plain data frame operation. Drop the geometry whenever you do not need it back. In the pipeline, county-joined-with-wells enters first, group-by on county-f-p defines one group per county code, and summarize creates mean-g-w-extraction by applying mean to the joined well values in each group. A county with no wells has only a missing extraction value, so its reported mean remains missing with the code exactly as shown. If you keep the geometry, s-f must union the repeated county geometries and the result stays spatial. If you uncomment st-drop-geometry, the geometry column disappears before grouping, the unnecessary union is avoided, and the result is a much faster ordinary table. You could replace the expression inside summarize to calculate a different county statistic without changing the spatial join.
Since we joined the two layers, we can now do calculations that were not possible before. Here, we will calculate the average groundwater extraction by county.
dplyr::summarize() takes a long time when it is applied to an sf object, because the geometry travels with the data and summarize() unions the geometries of each group as well as summarising the attributes.sf::st_drop_geometry() and it becomes an ordinary data.frame operation, which is far faster. Drop the geometry whenever you do not need it back.Of course, it is just as easy to get other types of statistics by simply modifying the summarize() part.
The third and messiest combination. Polygons against polygons is one-to-many like the last case, but with a complication that points do not have: overlap is a matter of degree. A county might lie almost entirely inside a watershed or barely clip its corner, and st_join treats those two cases identically. It records that they intersect and nothing more. That limitation is real, and this section builds towards it deliberately rather than hiding it, so keep it in mind as we work through the example. The final section of the deck is the fix.
For each of the observations (polygons) in the target data,
Two Iowa datasets that share no key at all. Nitrogen use is recorded by county; the hydrologic units, the watersheds, follow rivers and ignore county lines entirely. There is no column you could join these on. Plot them and you can see the problem: the two sets of boundaries cut across each other with no relationship whatsoever. The nitrogen figures are simulated, but the geography is real, and that mismatch between administrative boundaries and natural ones is among the most common reasons anyone needs a spatial join at all. Data is collected on the units that are convenient to collect on, not the units you want to analyse. Use both inner tabs to inspect what each row represents before joining. On Nitrogen use, printing i-a-nitrogen shows the county attributes and geometry, and the map passes that object to geom-s-f while mapping nitrogen-rate to fill, so colour displays the simulated pounds-per-acre values. On Hydrologic units, printing huc-i-a shows the HUC attributes and geometry, and its map draws those watershed polygons without an attribute fill. Flipping between the two tabs makes the incompatible boundaries, units of observation, and available attributes explicit. The Story tab then explains why those two layers need to meet.
Nitrogen use (lb/acre) by county in Iowa (Note: this is a fake dataset that is generated using R):
Hydrologic units that cover Iowa:
Here is the research question that motivates the whole section. You want to know whether nitrogen use affects water quality. Assume water quality is measured by watershed, while nitrogen use is recorded by county. This demonstration contains no water quality variable; it only attaches county nitrogen rates to HUC geometries. To run a regression with water quality data, you would need both variables on the same rows, and geography is the only thing that can put them there. Look at the overlaid map and think about what a sensible answer would even be for a watershed that spans eight counties. That question, what a sensible answer looks like when one watershed spans eight different counties, is exactly where this section ends up. The map puts that question on screen. The first geom-s-f draws the county layer and maps nitrogen-rate to fill; alpha at zero-point-six makes the fill translucent so boundaries from both systems can remain visible. The second geom-s-f adds the HUC layer with alpha equal to zero, leaving its outlines over the counties. Theme-void removes axes and grid lines, and the final theme call moves the nitrogen legend to the bottom where it does not cover the map. Follow those crossing outlines into Demonstration, where location supplies the missing connection.
You are interested in understanding the impact of nitrogen use for agricultural production on water quality.
huc_iaia_nitrogen)You would like to associate nitrogen use values with water quality values so that you can run statistical analysis on the impact of nitrogen use on water quality.
The join itself, one line as always, and then the row counts tell the story. The watershed with code 07060004 intersects eight counties, so it comes back as eight rows, each carrying the identical watershed geometry and one county’s nitrogen figure. Then the summarise gives you a mean per watershed. Notice that this treats all eight counties as equal, which is the assumption the next tab attacks. For now the mechanics are what matter: join, then group, then summarise, exactly the same three steps as the previous section. The assumption buried in that mean is what the next tab is about. In the st-join call, the huc-i-a polygons are the target because they are first, and the i-a-nitrogen polygons are the source because they are second. The result is saved as huc-joined-with-acres. Filtering where HUC-code equals the quoted code isolates the eight matching rows, and the surrounding parentheses print them for inspection. The summary pipeline then deliberately calls st-drop-geometry before grouping, because this output is a table of values and does not need an expensive geometry union. Group-by on HUC-code defines one watershed at a time, and mean of nitrogen-rate creates the simple average named average-nitrogen-rate. Move to But to see why that computationally valid mean is substantively weak.
Let’s join the two:
Here, for each of the HUC units from huc_ia, all the intersecting counties from ia_nitrogen are matched.
For example, HUC_CODE == "07060004" intersects eight counties.
All eight rows carry exactly the same geometry, the one representing the HUC unit with HUC_CODE == 07060004.
We can now find the average nitrogen use (lb/acre) by HUC unit:
The catch, and it is a serious one. The join told you that eight counties intersect this watershed and nothing whatsoever about how much of each. A county contributing ninety percent of the area and one clipping a corner both count once in that mean. Look at the map: the sizes of the overlaps are obviously wildly different. So a simple average of the eight is not a defensible number for a regression. The fix is area weighting, which needs a different function, and that is the final section of the deck. The important callout states the limitation precisely: the joined rows record the existence of an intersection, not its shape or size. In the plotting code, the filter uses the percent-in-percent operator to keep county FIPS codes found in huc-zero-seven-zero-six-zero-zero-zero-four, so only the eight matched counties are filled red. Their alpha of zero-point-five lets the HUC boundary remain visible, and the second layer draws that repeated HUC geometry on top. The picture shows why a binary match cannot supply weights. Follow the link to Cropping Join later in the deck, where the overlap geometry itself is returned and can be measured.
Important
Note that the resulting dataset does not tell you the nature of intersections. The only thing we know from huc_joined_with_acres is which counties the HUC units are intersecting with no matter how small or large the overlapping area are.
We simply take the average of the value of nitrogen_rate of the intersecting counties. But, this does not take into account the degree of the overlaps between the intersecting counties and the HUC unit. Later, we will talk about how to find area-weighted average of nitrogen_rate this section.
Intersects is the default relation but not the only one. The join argument takes a function describing whichever relation you want. Additional arguments that relation needs, such as a distance, can be supplied directly to st_join, which forwards them to the relation. For example, st_join with join set to st_is_within_distance and dist set to five is valid. Wrapping the relation and its distance in a little function is optional. The next few slides put the distance relation to work on a problem that intersects cannot solve at all.
Spatial join with st_join() uses st_intersects() as the default topological relationship for joining. You can pick a different one with the join argument.
Syntax
where st_* determines the topological relationship between sf_1 and sf_2, and ... supplies any additional arguments that function needs.
Extra arguments are forwarded to the relation
st_join(sf_1, sf_2, join = st_is_within_distance, dist = 5) is valid because st_join() forwards dist to the join function. Writing join = \(x, y) st_is_within_distance(x, y, dist = 5) is also valid, but optional.
Two point layers from an on-farm experiment: soybean yield measured by the combine, and seed rate recorded by the planter. Both are dense clouds of GPS readings over the same field. The reason this needs a distance relation rather than intersects is worth pausing on. These are two machines making separate passes over the same ground, so no yield point ever sits exactly on a seed rate point. Intersects would match nothing whatsoever, and you would get a dataset of missing values with no error to tell you why. Proximity is the only relation that can work here. The two data calls load soy-yield and as-applied-s-rate from the course data package into the working session. Both arrive as s-f point objects in the same projected coordinate reference system, whose distance unit is metres. That shared projected C-R-S is why the distance threshold used later can be interpreted directly in metres rather than degrees. The next tab supplies the research objective and lets you see how the two point clouds relate.
We use soy bean yield (points) data and as-applied seed rate (points) data.
The economic question is whether planting more seed raises yield, and answering it means having both numbers on the same row. Look at the plot of red and blue points and you can see they interleave without ever coinciding. So the rule becomes: for each yield point, take the seed rate points within ten metres. That threshold is a judgement call about how far apart two readings can be and still describe the same patch of ground, and it is the sort of choice you should be prepared to defend. The empty ggplot call starts a layered map. The first geom-s-f draws soy-yield as small red points, with size set to zero-point-four so the dense readings do not obscure one another. The second uses the same size for as-applied-s-rate but colours those points blue, making the two machines’ tracks distinguishable. Theme-void removes the axes and leaves the spatial pattern. The plot justifies the proximity rule visually; the Demonstration tab turns that rule into the join argument.
You have run an on-farm field experiment to understand the impact of seed rate on soybean yield. Both are available as points data.
You want to merge them together based on their proximity so that you can run statistical analysis. Specifically, for each of the yield points, we would like to link the seed rate points that are within 10-meter from the yield point.
The join with a distance relation, written as a small anonymous function with the ten metre threshold inside it. Then the summarise, because this is one-to-many again and a yield point may match several seed rate points. Note the commented-out line dropping the geometry, which is the same performance point from earlier and applies just as much here. The result is one average seed rate per yield point, which is the dataset you would actually run a regression on. Everything before this was getting the two machines’ readings onto the same rows. Read the call carefully. Soy-yield is first, so each yield observation is a target; as-applied-s-rate is second, so seed attributes are the source. The anonymous function receives the two geometry sets as x and y and returns the matches from st-is-within-distance; dist equals ten in the layers’ metre-based C-R-S. The joined one-to-many rows are saved as soy-seed. In the next pipeline, group-by on yield-i-d gathers all matches for each original yield reading, and mean of seed-rate names their average average-seed-rate. An unmatched yield point still produces a missing average. Uncommenting st-drop-geometry is appropriate when the regression table no longer needs point geometry. Move to Inspect and verify before trusting the threshold.
Let’s join using st_join():
We can now summarize the joined data like below:
The verification, and this is the tab worth spending time on. The first yield point matched nothing, so its seed rate is missing, which happens whenever nothing lies within the threshold. The second matched two points. Rather than take that on trust, the code draws it: the yield point in red, a ten metre circle around it, and the nearby seed rate points labelled with their identifiers. Count what falls inside the circle. This is the same verify-by-picture habit from the first join, and it is how you catch a wrong threshold or a wrong CRS. Printing soy-seed exposes the repeated matches and missing values directly. For the picture, the spatial subset uses a 20-metre buffer around row two of soy-yield only to collect enough nearby seed points to label; it does not change the join’s 10-metre rule. The plot draws that second yield point in red, then draws the actual 10-metre buffer as an unfilled blue boundary. Geom-s-f-text places each nearby point’s seed-i-d at its location with size six, and theme-void removes distractions. IDs one and five-five-eight fall inside the blue circle, so the visual check agrees with the two rows produced by the join.
According to the join, the 1st yield point from soy_yield did not have any seed rate data points from as_applied_s_rate within its 10 meter radius, so NA in seed_rate.
The second yield point from soy_yield is matched with two seed rate points from as_applied_s_rate: seed_id = 1 and seed_id = 558. Are they indeed less than 10 metres from the second yield point?
aggregate()A shortcut. Every one-to-many example so far took two steps, join and then summarise, and aggregate does both at once. Watch the argument order, because it is the opposite of st_join and that is a genuine trap: here the thing being aggregated comes first and the thing you are grouping by comes second. The function you pass as FUN is applied to every column of the source layer, which is convenient and occasionally silly, as the demonstration on the next tab shows. Intersects is the default relation here too, and as with st_join you can change it.
In the example of finding average groundwater use by county (go here to remind yourself of this example), we took a two-step procedure to do so.
st_join()dplyr::summarize() to the joined objectHowever, this can actually be done in one step using aggregate(), in which you specify how you want to aggregate with the FUN option:
Syntax
Here, for each of the rows of the second sf (here, polygons_sf), all the intersecting points in points_sf are found and then the average of all the columns of points_sf are calculated. Yes, st_intersects() is the default topological relationship just like spatial sub-setting with sf1[sf2, ] and st_join().
The one-liner, and then the reason to be careful with it. Look at the output: the well identifier has been averaged along with the extraction figures, giving a meaningless number in a column that looks perfectly respectable. That is what applying one function to every column gets you. The fix is on the same slide, selecting the column you actually want before aggregating. The general lesson is worth keeping: an operation that quietly does something to every column will do it to the columns you were not thinking about. In the first call, wells-ne supplies the values to aggregate, ne-counties supplies one output geometry per county, and FUN equals mean tells aggregate how to combine every matched source column. The second call wraps the first argument in d-plyr select, keeping g-w-extraction as the only attribute to average. Because this is an s-f object, its geometry is retained even though it is not named in select, so aggregate can still find which wells intersect each county. The resulting county layer now contains a meaningful mean extraction column without the meaningless mean well identifier.
Note that wellid was also averaged by county. We could just do this:
The same shortcut applied to the soybean example, which shows that aggregate is not limited to intersects — the join argument works here exactly as it did with st_join. Note the warning in the bullets about argument order, because the code looks almost identical to the earlier version while meaning something different. Also be aware that this cell overwrites the soy_seed object you built a few slides ago with a different result, so if you go back and re-read that tab, you will not be looking at what it made. Here as-applied-s-rate comes first because its columns are being averaged, while soy-yield comes second because each yield-point geometry defines an output group. FUN equals mean applies that average, and the anonymous join function passes x, y, and a 10-metre dist to st-is-within-distance. The outer parentheses print the newly assigned result. It has one geometry for each yield point, including missing source values where nothing matched. It also averages seed-i-d as well as seed-rate, so the warning from the preceding tab still applies: select only the source attributes whose means have a sensible interpretation.
aggregate() is a fairly general procedure of spatial joining and summarization, and you can use it for many other cases including our example of soybean yield and seed rate (go here for the example).
Here is the code:
This looks almost identical with the code to spatial join the two sf layers. However, the order of the sf objects is flipped.
You are aggregating the columns of first sf for each of the geometry in the second sf (here yield point).
Back to the problem the polygons-against-polygons section left open. A simple average across eight counties ignores how much of each one actually falls inside the watershed, and we want a weighted one. The insight is that if you can get the geometry of each overlapping piece, you can measure its area, and area is exactly the weight you need. Everything left in this deck is about producing those pieces of geometry, and it turns out that a single function does it. The weighting itself is then just arithmetic you already know how to do.
In the example of finding average nitrogen use by HUC unit using st_join(), we had a problem of not knowing how much of each of the intersecting counties shares the same area with the HUC unit.
If we can get the geometry of the intersecting part of the HUC unit and the county, then we can calculate its area, which in turn allows us to find area-weighted averages of joined attributes.
The function is st_intersection, and the name is confusingly close to st_intersects, so hold the difference clearly. Intersects, ending in s, asks a yes-or-no question and returns indices. Intersection, ending in tion, actually cuts the geometry and hands you back the overlapping pieces themselves, with the non-overlapping parts removed. It also carries the attributes across from both layers. That combination — real geometry you can measure, plus the attributes of both layers on the same row — is precisely what makes area weighting possible. The next two tabs show it on data small enough to check by eye.
For these purposes, we can use sf::st_intersection().
While st_intersects() returns the indices of intersecting objects, st_intersection() returns intersecting spatial objects with the non-intersecting parts of the sf objects cut out.
Moreover, attribute values of the source sf will be merged to its intersecting sfg in the target sf.
Back to the small made-up data from last deck, because with three polygons and two lines you can check the answer by eye. Intersecting the lines with the polygons gives you the line segments that fall inside each polygon, and every segment is one row labelled with which line and which polygon produced it. Compare the two figures side by side: on the left the full lines running out past the polygon edges, on the right only the parts inside, cut cleanly at the boundary. That cutting is the whole difference from subsetting. The setup code constructs each geometry from coordinates. St-point makes the three individual points. Each st-linestring turns a matrix of endpoint coordinates into a line. Each polygon is a closed coordinate ring whose final coordinate repeats the first, wrapped in a list for st-polygon. For each collection, list gathers the geometries, st-s-f-c makes a simple-feature geometry column, st-as-s-f makes an s-f object, and mutate adds readable point, line, or polygon names. The points are part of the reusable toy setup even though this particular comparison uses only lines and polygons. The active operation is st-intersection of lines and polygons, so it returns every portion of a line that lies within a polygon and carries both layers’ name attributes. The pipe adds int-name by using paste-zero to join the line name, a hyphen, and the polygon name; the surrounding parentheses both save and print intersections-l-p. In the left plot, polygon name controls fill with alpha zero-point-three, line name controls colour, and the two discrete scales give the legends the titles Polygons and Lines. In the right plot, int-name controls the clipped segments’ colours, the translucent polygons provide context, and theme-void keeps both maps focused on geometry. The next tab repeats the operation with polygon targets and sources.
The following code gets the intersection of the lines and the polygons.
Here is how the lines and polygons look like:
Here is how the intersections and polygons look like.
The same operation with polygons, which is the case we actually need. Intersecting polygons one and three with polygon two returns one row for each overlap, with attributes from both parents. Polygon one is clipped into a new geometry, while polygon three remains unchanged because it lies wholly inside polygon two. Look at the figure: the original polygons are filled and faint, and the intersections are outlined in heavy colour, sitting exactly where the shapes overlap. Crucially, those outlined pieces have areas you can measure. That is what we came for, and the next tab applies it to the watersheds. The c of one and three row subscript selects polygons one and three as the first input, while the row-two subscript selects polygon two as the second. St-intersection cuts each first-input geometry to its overlap with the second and combines their attributes. Because both inputs contain a column called polygon-name, s-f keeps the first name unchanged and gives the second one the suffix dot-one. Paste-zero combines those two columns with a hyphen to form int-name, and the outer parentheses save and print intersections-p-p. The verification plot first maps every original polygon name to fill and uses alpha zero-point-three so overlaps remain visible. A second geom-s-f draws the returned pieces, mapping int-name to outline colour, increasing line-width to one-point-five, and setting alpha to zero so the new fill does not hide the originals. Theme-void removes the coordinate furniture. The coloured outlines are therefore the measurable geometries returned by st-intersection, not merely a highlight applied to the originals.
The following code gets the intersection of polygon 1 and polygon 3 with polygon 2. Each instance of the intersections of polygons 1 and 3 against polygon 2 becomes an observation (polygon 1-polygon 2 and polygon 3-polygon 2).
Just like the lines-polygons case, the non-intersecting part of polygons 1 and 3 are cut out and do not remain in the returned sf.
And here it is applied to the real problem. Intersecting the watersheds with the counties gives one row per watershed-county pair, each with the geometry of that specific overlap. Take the area of each piece, multiply each county’s nitrogen rate by its area, add those up and divide by the total area, and you have a properly weighted average in which a county contributing ninety percent of the watershed counts ninety times as much as one contributing one percent. The map at the end shows the pieces one watershed is made of. In the first pipeline, st-intersection with huc-i-a first and i-a-nitrogen second cuts each HUC by every county it overlaps and carries the HUC and nitrogen attributes onto each returned piece. Mutate then creates huc-county by concatenating HUC-code, a hyphen, and county-f-p, giving each pair a useful label. The surrounding parentheses save the object as HUC-intersections and print it. Unlike st-join, its repeated HUC rows no longer carry identical whole-watershed geometries. Each row now contains only one HUC-county overlap. The final pipeline filters to the quoted HUC code zero-seven-zero-six-zero-zero-zero-four, passes those eight pieces into ggplot, and maps huc-county to fill so each county contribution has a different colour. Theme-void removes the axes. This tab prepares and displays the pieces but does not yet calculate their areas. The Exercise tab gives you the exact st-area, multiplication, sum, and division pattern needed to finish that calculation.
Let’s now get back to the example of HUC units and county-level nitrogen use data. We would like to find area-weighted average of nitrogen use instead of the simple average.
Using st_intersection(), for each of the HUC polygons, we can find the intersecting counties, and then divide it into parts based on the boundary of the intersecting polygons.
The key difference from the st_join() example is that each observation of the returned data is a unique HUC-county intersection.
The figure below maps every intersection of the HUC unit HUC_CODE == 07060004 with the eight counties it overlaps.
Your turn, on a deliberately simple layout: four equal two-by-two squares in a grid with values one to four, and another two-by-two square sitting over the middle of all four. Because it is centred, each overlap is one unit square, so you can work out the weighted answer in your head before writing any code. That is the point of the setup. Do problem one to get the intersection, then problem two to weight by area, and check your result against the number you predicted. If they disagree, the code is wrong, not your arithmetic. Start in Data Preparation. Each st-polygon receives a list containing a closed matrix of corner coordinates, and st-s-f-c combines the four adjacent squares into one geometry column. St-s-f with geometry equal to dot turns that column into an s-f object, and mutate with value equal to one through four assigns the four values in row order. The second construction uses the same steps for the single red square from coordinates one-one to three-three. Printing both objects lets you inspect attributes and geometry before plotting. In the map, factor of value makes fill a four-category legend rather than a continuous scale, alpha zero-point-five lets the layers show through one another, and theme-void leaves only the geometry. For Problem 1, fill in the editable st-intersection call with the two objects. The folded answer uses polygon-two first and polygon-one second, then saves the four one-by-one overlap pieces as intersection. For Problem 2, st-area of geometry measures each returned piece. As-numeric removes the units wrapper so those areas can be used in the arithmetic shown. Mutate stores the result as area, and summarize computes the weighted mean as the sum of value times area divided by the sum of area. Here every area is one, so the expression returns two-point-five, the ordinary mean of one, two, three, and four. In a real watershed, unequal overlap areas are exactly what make this formula different from a simple mean.
We use polygon_1 and polygon_2 in this exercise. Inspect them to familiarize yourself with them:
Here is what they look like:
Find the intersection of the two polygons using st_intersection().
Answer