R as GIS: Interaction of Vector Datasets I

Before you start


Learning objectives

The objectives of this chapter is to learn spatial interactive operations that involves two sf objects. Specifically,

  • understand topological relations
  • subsetting an sf object based on another sf object


Tips to make the most of the lecture notes


Interactive navigation tools

  • 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


Running and writing codes

  • The box area with a hint of blue as the background color is where you can write code (hereafter referred to as the “code area”).
  • Hit the “Run Code” button to execute all the code inside the code area.
  • You can evaluate (run) code selectively by highlighting the parts you want to run and hitting Command + Enter for Mac (Ctrl + Enter for Windows).
  • If you want to run the codes on your computer, you can first click on the icon with two sheets of paper stacked on top of each other (top right corner of the code chunk), which copies the code in the code area. You can then paste it onto your computer.
  • You can click on the reload button (top right corner of the code chunk, left to the copy button) to revert back to the original code.

Topological Relations

Definition

Topological relations refer to the way multiple spatial objects are spatially related to one another.


Goals

You can identify various types of spatial relations using the sf package

  • Our main focus is on the intersections of spatial objects, which can be found using st_intersects().

  • We also briefly cover st_is_within_distance() and st_nearest_feature()

You can run ?geos_binary_pred to find all the topological relations you can use:

We use three sf objects: points, lines, and polygons. Inspect each of them.



If you are interested in codes to create them see below.

Code
#--- create points ---#
point_1 <- sf::st_point(c(2, 2))
point_2 <- sf::st_point(c(1, 1))
point_3 <- sf::st_point(c(1, 3))

#--- combine the points to make a single  sf of points ---#
points <- list(point_1, point_2, point_3) %>% 
  sf::st_sfc() %>% 
  sf::st_as_sf() %>% 
  mutate(point_name = c("point 1", "point 2", "point 3"))

#--- create lines ---#
line_1 <- sf::st_linestring(rbind(c(0, 0), c(2.5, 0.5)))
line_2 <- sf::st_linestring(rbind(c(1.5, 0.5), c(2.5, 2)))

#--- combine the points to make a single  sf of points ---#
lines <- list(line_1, line_2) %>% 
  sf::st_sfc() %>% 
  sf::st_as_sf() %>% 
  mutate(line_name = c("line 1", "line 2"))

#--- create polygons ---#
polygon_1 <- sf::st_polygon(list(
  rbind(c(0, 0), c(2, 0), c(2, 2), c(0, 2), c(0, 0)) 
))

polygon_2 <- sf::st_polygon(list(
  rbind(c(0.5, 1.5), c(0.5, 3.5), c(2.5, 3.5), c(2.5, 1.5), c(0.5, 1.5)) 
))

polygon_3 <- sf::st_polygon(list(
  rbind(c(0.5, 2.5), c(0.5, 3.2), c(2.3, 3.2), c(2, 2), c(0.5, 2.5)) 
))

#--- combine the polygons to make an sf of polygons ---#
polygons <- list(polygon_1, polygon_2, polygon_3) %>% 
  sf::st_sfc() %>% 
  sf::st_as_sf() %>% 
  mutate(polygon_name = c("polygon 1", "polygon 2", "polygon 3"))

st_intersects() checks which of sfgs in an sf geographically intersect with which of sfgs in another sf.



  • The output is a list of which polygon(s) each of the points intersect with.

  • The numbers 1, 2, and 3 in the first row mean that 1st (polygon 1), 2nd (polygon 2), and 3rd (polygon 3) objects of the polygons intersect with the first point (point 1) of the points object.

  • The fact that point 1 is considered to be intersecting with polygon 2 means that the area inside the border is considered a part of the polygon (of course).


If you would like the results of st_intersects() in a matrix form with boolean values filling the matrix, you can add sparse = FALSE option.


The output is a list of which polygon(s) each of the lines intersect with.


For polygons vs polygons interaction, st_intersects() identifies any polygons that either touches (even at a point like polygons 1 and 3) or share some area.

sf::st_is_within_distance() function identifies whether any of sfgs in sf_2 is within the specified distance from each of the sfgs in sf_1.

st_is_within_distance(sf_1, sf_2)

Create two sets of points and then inspect each of them.


Here is the visualization of the two sets of points we just created.

We want to know which of the blue points (points_set_2) are located within 0.2 from each of the red points (points_set_1).

The following figure gives us the answer visually.



Confirm that your visual inspection results are consistent with the outcome of the following code using st_nearest_feature() function.

sf::st_nearest_feature() identifies which sfgs in sf_2 is closest in distance to each of sf_1.

st_nearest_feature(sf_1, sf_2)

Confirm that your visual inspection results are consistent with the outcome of the above code using st_is_within_distance() function.

Run the following codes to get mower_sensor and fairway_grid data.


Visualize the datasets.

Use sf::st_intersects() to find out which of the points in mower_sensor_sf are inside of any of the polygons in fairway_grid.


Since there are so many points in mower_sensor_sf, you won’t really see which ones are inside of any of the polygons in fairway_grid.

That is okay for now. We will later learn how to filter sf objects spatially.


Answer codes

Code
st_intersects(mower_sensor_sf, fairway_grid)

Spatial Subsetting

Spatial subsetting refers to operations that narrow down the geographic scope of a spatial object based on another spatial object.

Here are the datasets we will use here. Inspect them to familiarize yourself with the datasets.

Select only the counties that intersect with the HPA boundary.

When subsetting a data.frame by specifying the row numbers you would like to select, you can do

#--- NOT RUN ---#
data.frame[vector of row numbers, ]


Spatial subsetting of sf objects works in a similar syntax:

#--- NOT RUN ---#
sf_1[sf_2, ]


where you are subsetting sf_1 based on sf_2. Instead of row numbers, you provide another sf object in place.

The following code spatially subsets Nebraska counties based on the HPA boundary.


  • You can see that only the counties that intersect with the HPA boundary remained.

  • This is because when you use the above syntax of sf_1[sf_2, ], the default underlying topological relation is st_intersects().

  • So, if an object in sf_1 intersects with any of the objects in sf_2 even slightly, then it will remain after subsetting.

Sometimes, you just want to flag whether two spatial objects intersect or not, instead of dropping non-overlapping observations like we saw with sf_1[sf_2, ] syntax. In that case, you can get a list of the IDs and then assign 1 (or TRUE) if in the list, 0 (or FALSE) otherwise.



Get the list of countyfp (ID) of the intersected counties:


Assign 1 or 0 to a new variable called in_hpa based on the list.

You can specify the topological relation as in

#--- NOT RUN ---#
sf_1[sf_2, op = topological_relation_type] 


For example, if you only want counties that are completely within the HPA boundary, you can do the following:


Check visually:

Select only the wells that intersect with (or equivalently inside) the HPA boundary.

We can select only the wells that reside within the HPA boundary using the same syntax as the polygon-polygon example.


As you can see in the figure below, only the wells that are intersects the HPA remained because the default topological relation is st_intersects() (here, you will get the same results even if you use op = st_within.).

Get the list of wellid (ID) of the intersected wells:


Assign 1 or 0 to a new variable called in_hpa based on the list.

Select only railroads that intersect with the Lancaster county in Nebraska.


Just like we did in the two previous examples:


Successful?

Get the list of LINEARID (ID) of the intersected railroads:


Assign 1 or 0 to a new variable called in_hpa based on the list.

Spatial Cropping

  • We can use st_crop() to crop spatial objects to a spatial bounding box (extent) of a spatial object.

  • The bounding box of an sf is a rectangle represented by the minimum and maximum of x and y that encompass/contain all the spatial objects in the sf.

  • You can use st_bbox() to find the bounding box of an sf object.

Let’s get the bounding box of the High-Plains aquifer using st_bbox().


Check its class:


You can convert a bbox to sfc by applying st_as_sfc() to the bbox object (you cannot use a bbox for mapping and other interactive spatial operations).


The bounding box looks like this (red rectangle):

`

Now, let’s crop Nebraska counties to the bounding box of the High-Plains aquifer boundary.


Note that, you do not need to do the following (they would produce the same outcome):


Note that st_crop() will chop off the parts that are not intersecting.