Spatial Queries: Find Yosemite POIs in the Upeer Merced Subbasin

In this Notebook we’ll use geoprocessing functions from sf to identify Yosemite Points-of-Interest that fall within the Upper Merced Subbasin.

Setup

Load the packages we’ll need and set tmap mode to ‘plot’:

library(sf)
library(tmap)
tmap_mode("plot")

Load dplyr and set name conflict preferences:

library(dplyr)

## Load the conflicted package
library(conflicted)

# Set conflict preferences
conflict_prefer("filter", "dplyr", quiet = TRUE)
conflict_prefer("count", "dplyr", quiet = TRUE)
conflict_prefer("select", "dplyr", quiet = TRUE)
conflict_prefer("arrange", "dplyr", quiet = TRUE)

Practice Querying with Sample Data

Import Practice Data

First we import some practice data:

circles_sf <- st_read("./data/test_circles.geojson")

Reading layer `test_circles' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data\test_circles.geojson' using driver `GeoJSON'
Simple feature collection with 3 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -0.5 ymin: -0.5 xmax: 4 ymax: 3
Projected CRS: WGS 84 / UTM zone 11N

circles_sf

Simple feature collection with 3 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -0.5 ymin: -0.5 xmax: 4 ymax: 3
Projected CRS: WGS 84 / UTM zone 11N
  circle_id                       geometry
1         A POLYGON ((1.5 0.5, 1.499391...
2         B POLYGON ((4 0.5, 3.999391 0...
3         C POLYGON ((4 2, 3.999391 2.0...

pts_sf <- st_read("./data/test_pts.geojson")

Reading layer `test_pts' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data\test_pts.geojson' using driver `GeoJSON'
Simple feature collection with 120 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -0.47835 ymin: -0.4040069 xmax: 3.933032 ymax: 2.985598
Projected CRS: WGS 84 / UTM zone 11N

pts_sf

Simple feature collection with 120 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -0.47835 ymin: -0.4040069 xmax: 3.933032 ymax: 2.985598
Projected CRS: WGS 84 / UTM zone 11N
First 10 features:
   pt_id                     geometry
1      1 POINT (-0.4430744 0.9685033)
2      2  POINT (0.3676537 0.2035441)
3      3   POINT (1.535087 0.9823107)
4      4   POINT (0.8340309 2.645155)
5      5     POINT (2.901864 1.73253)
6      6   POINT (3.663181 0.7940655)
7      7   POINT (0.8618392 2.640156)
8      8    POINT (2.832866 2.870026)
9      9   POINT (3.046212 0.8546275)
10    10  POINT (-0.349015 0.7595088)

Plot the points on top of the circles:

tm_shape(circles_sf) +
  tm_borders(col = palette()[2:4] ) +
  tm_text("circle_id") +
tm_shape(pts_sf) +
  tm_dots(col = "dimgray") +
tm_grid(labels.show = TRUE, lines = FALSE)

Identify the points in Circle A

Next we identify the points in circle A using a spatial predicate function (st_intersects). We could also copy the points in Cirlce A with st_intersection(), but there are times when you don’t need or want to make copies of the data.

circle_a_sf <- circles_sf %>% filter(circle_id == "A")

pt_in_circle_a_yn_mat <- pts_sf %>% 
  st_intersects(circle_a_sf,
                sparse = FALSE)
  
head(pt_in_circle_a_yn_mat)

      [,1]
[1,] FALSE
[2,]  TRUE
[3,] FALSE
[4,] FALSE
[5,] FALSE
[6,] FALSE

Since we have a column of TRUE/FALSE values, we can subset those features using the dplyr filter function:

## Copy the points in circle A to a new object. Note in the filter expression
## we use square bracket notation to pull out the first column of the matrix 

a_pts_sf <- pts_sf %>% 
  filter(pt_in_circle_a_yn_mat[,1])

## Plot to verify
tm_shape(circles_sf) +
  tm_borders(col = palette()[2:4] ) +
  tm_text("circle_id") +
tm_shape(a_pts_sf) +
  tm_dots(col = "red", size = 0.1) +
tm_grid(labels.show = TRUE, lines = FALSE)

CHALLENGE: How many points in cirlce B?

Answer

circle_b_sf <- circles_sf %>%
  filter(circle_id == "B")

pts_sf %>% 
  st_intersects(circle_b_sf, sparse = FALSE) %>% 
  sum()

[1] 26

CHALLENGE: Plot the points that fall within Circle B and Circle C

Answer

pts_circle_bc_sf <- pts_sf %>% 
  st_intersection(circles_sf %>% filter(circle_id == "B")) %>% 
  st_intersection(circles_sf %>% filter(circle_id == "C"))

attribute variables are assumed to be spatially constant throughout all geometriesattribute variables are assumed to be spatially constant throughout all geometries

tm_shape(circles_sf) +
  tm_borders(col = palette()[2:4] ) +
  tm_text("circle_id") +
tm_shape(pts_circle_bc_sf) +
  tm_dots(col = "red", size = 0.1) +
tm_grid(labels.show = TRUE, lines = FALSE)

CHALLENGE: How many points don’t fall in any circle?

This is a little more challenging. One approach we could take is to test for intersection for all points and all circles. This will return a matrix of TRUE / FALSE values, with one row for each point and one column for each circle. If a row has three FALSE values, it means that point doesn’t intersect any circle.

Answer

## Approach one: Test intersection for all circles
intersect_yn_mat <- pts_sf %>% 
  st_intersects(circles_sf, sparse = FALSE)
head(intersect_yn_mat)

      [,1]  [,2]  [,3]
[1,] FALSE FALSE FALSE
[2,]  TRUE FALSE FALSE
[3,] FALSE FALSE FALSE
[4,] FALSE FALSE FALSE
[5,] FALSE FALSE  TRUE
[6,] FALSE  TRUE FALSE

## Sum up the rows. If a the point doesn't intersect any circle, the sum of the row will be 0
outside_all_circles <- rowSums(intersect_yn_mat)

## Count the number of rows where the total is 0. These are the points that don't intersect any circle.
sum(outside_all_circles == 0)

[1] 51

## Plot to make sure
pts_outside_circles_sf <- pts_sf %>% 
  filter(outside_all_circles == 0)

tm_shape(circles_sf) +
  tm_borders(col = palette()[2:4] ) +
  tm_text("circle_id") +
tm_shape(pts_outside_circles_sf) +
  tm_dots(col = "dimgray", size = 0.1) +
tm_grid(labels.show = TRUE, lines = FALSE)

A second approach would be to union the circles and then take the difference:

circle_union_sf <- circles_sf %>% st_union(by_feature = FALSE)

pts_outside_circle_union_sf <- pts_sf %>% 
  st_difference(circle_union_sf)

attribute variables are assumed to be spatially constant throughout all geometries

tm_shape(circle_union_sf) +
  tm_fill() +
tm_shape(pts_outside_circle_union_sf) +
  tm_dots(col = "orange", size = 0.15) +
tm_grid(labels.show = TRUE, lines = FALSE)

CHALLENGE: Plot the points that lie within 0.25 map units of a circle, but are not contained within the circle

One approach we can use to solve this problem is to feed the point layer into st_is_within_distance() (which returns a matrix), then feed that matrix into rowSums() (which adds up the rows converting TRUE to 1 and FALSE to 0), and then feed the row sums into a condition check (row sum > 0). This will give us a logical vector (TRUE / FALSE) values that indicate whether each point is close to a circle.

Unfortunately st_is_within_distance() also returns TRUE when a point is within the circle (because distance = 0). Hence we’ll need to remove the points within the circles from our set. We can do this by creating a logical vector for intersects, and then write a logical expression that uses the ‘and’ and ‘not’ operators to find the points we want.

Answer

library(magrittr)

## Test whether points are near a circle. Note st_is_within_distance() also returns
## TRUE for points that are within the circle (because the distance = 0)
near_or_within_any_circle_yn <- pts_sf %>% 
  st_is_within_distance(circles_sf, dist = 0.25, sparse = FALSE) %>% 
  rowSums() %>% 
  is_greater_than(0)    ## you could also use {. > 0}

## Test the points that are within circles (so we can remove them)
within_any_circle_yn <- pts_sf %>% 
  st_intersects(circles_sf, sparse = FALSE) %>% 
  rowSums() %>% 
  is_greater_than(0)    ## you could also use {. > 0}

## Use the logical 'and' operator (&) and 'not' operator (!) to get those
## rows near but not within circles
near_but_not_within_yn <- near_or_within_any_circle_yn & !within_any_circle_yn

## Count the TRUEs
sum(near_but_not_within_yn)

[1] 21

## Plot to make sure
tm_shape(circles_sf) +
  tm_borders(col = palette()[2:4] ) +
  tm_text("circle_id") +
tm_shape(pts_sf %>% filter(near_but_not_within_yn)) +
  tm_dots(col = "brown", size = 0.1) +
tm_grid(labels.show = TRUE, lines = FALSE)

Another approach we can take is to create a ‘donut ring’ around the union of the threes circles, and then use st_intersecs() or st_intersection() to plot the points we want. You can make a donut ring with a combination of st_buffer() and st_difference().

## Merge the circles into a single multipolygon feature
circles_merged_sf <- circles_sf %>% st_union()

## Create the donut ring:
circles_ring_sf <- circles_merged_sf %>% 
  st_buffer(dist = 0.25) %>% 
  st_difference(circles_merged_sf)
  
## Plot to make sure
tm_shape(circles_ring_sf) +
  tm_polygons()

## Find and plot the points that fall within the donut ring
pts_ring_sf <- pts_sf %>% st_intersection(circles_ring_sf)

attribute variables are assumed to be spatially constant throughout all geometries

## Count them
nrow(pts_ring_sf)

[1] 21

## Plot to verify
tm_shape(circles_sf) +
  tm_borders(col = palette()[2:4] ) +
  tm_text("circle_id") +
tm_shape(pts_ring_sf) +
  tm_dots(col = "purple", size = 0.1) +
tm_grid(labels.show = TRUE, lines = FALSE)

Plot the Points of Interest that Fall within the Upper Merced Subbasin

Next, we’ll apply what we learned to find the Yosemite Points-of-Interest that fall within the Upper Merced HUB-8 Subbasin.

Import the Watersheds

Start by importing the planning watershed units from calw221:

## Import the planning watersheds
gpkg_watershd_fn <- "./data/yose_watersheds.gpkg"
yose_watersheds_sf <- st_read(gpkg_watershd_fn, layer="calw221") 

Reading layer `calw221' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data\yose_watersheds.gpkg' using driver `GPKG'
Simple feature collection with 127 features and 10 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 1383.82 ymin: -61442.93 xmax: 81596.71 ymax: 26405.66
Projected CRS: NAD83 / California Albers

## View attribute table
yose_watersheds_sf %>% st_drop_geometry() %>% slice(1:6)

## Plot results
tmap_mode("plot")

tmap mode set to plotting

tm_shape(yose_watersheds_sf) + 
  tm_polygons("MAP_COLORS", palette = "Pastel1") 

Note the keyword MAP_COLORS tells tmap to select colors at random such that adjacent polygons have different colors.

Lump the Planning Watersheds into HUC-8 Subbasins

Next we’ll group the little planning watersheds into bigger “HUC-8” subbasins. This is easy because there is a column for the HUC 8 id number (HUC_8) and name (HUC_8_NAME).

yose_huc8_sf <- yose_watersheds_sf %>% 
  group_by(HUC_8) %>% 
  summarise(HUC_8_NAME = first(HUC_8_NAME), num_pws = n())

yose_huc8_sf

Simple feature collection with 6 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 1383.82 ymin: -61442.93 xmax: 81596.71 ymax: 26405.66
Projected CRS: NAD83 / California Albers

A convenient feature of group_by() is that when applied to a simple feature data frame it will also spatially aggregate (i.e., union) the features based on common values in the grouping column. Plot to verify:

epsg_utm11n_nad83 <- 26911
yose_bnd_utm <- st_read(dsn="./data", layer="yose_boundary") %>% 
  st_transform(epsg_utm11n_nad83)

Reading layer `yose_boundary' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 11 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -119.8864 ymin: 37.4947 xmax: -119.1964 ymax: 38.18515
Geodetic CRS:  North_American_Datum_1983

tm_shape(yose_huc8_sf) + 
  tm_polygons("MAP_COLORS", palette = "Pastel1") +
  tm_text("HUC_8_NAME", size = 0.7) +
tm_shape(yose_bnd_utm) +
  tm_borders(col = "red", lwd = 2)

Extract Upper Merced HUC-8 Subbasin

Next, pull out just the Upper Merced subbasin and save it as a separate object:

## Filter out just the Upper Merced Subbasin
merced_huc8_sf <- yose_huc8_sf %>% 
  filter(HUC_8_NAME == "UPPER_MERCED")
merced_huc8_sf

Simple feature collection with 1 feature and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2242.555 ymin: -61305.46 xmax: 65367.14 ymax: -12528.6
Projected CRS: NAD83 / California Albers

Import the Points-of-Interest

Import the Yosemite POIs:

## Import points of interest
yose_poi_utm <- st_read(dsn="./data", layer="yose_poi") %>% 
  select(OBJECTID, POINAME, POILABEL, POITYPE)

Reading layer `yose_poi' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data' using driver `ESRI Shapefile'
Simple feature collection with 2720 features and 30 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 246416.2 ymin: 4153717 xmax: 301510.7 ymax: 4208419
Projected CRS: NAD83 / UTM zone 11N

Identify Intersecting Points-of-Interest

Find out which POIs intersect the Upper Merced subbasin with st_intersects():

try(merced_poi <- yose_poi_utm %>% st_intersects(merced_watershed))

Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared,  : 
  object 'merced_watershed' not found

Oh no - ERROR message! Spatial querying requires features to be in the same CRS!

To fix this, we can project the Merced HUC-8 layer (which is in CA Albers) to match the POIs (which are UTM):

merced_huc8_utm_sf <- merced_huc8_sf %>% 
  st_transform(st_crs(yose_poi_utm))

Try the intersection test again:

yose_poi_merced_mat <- yose_poi_utm %>% st_intersects(merced_huc8_utm_sf, sparse=FALSE)
head(yose_poi_merced_mat)

     [,1]
[1,] TRUE
[2,] TRUE
[3,] TRUE
[4,] TRUE
[5,] TRUE
[6,] TRUE

CHALLENGE: How many points-of-interest fall within the Upper Merced subbasin?

Hint 1: This is equivalent to asking how many TRUE values there are in the first column of yose_poi_merced_mat.

Hint 2: To get the first column of a matrix x, use x[ , 1].

Answer

## How many points of interest fall in the Upper Merced subbasin?

## Option 1: Sum up the values in the column 1 of the matrix
sum(yose_poi_merced_mat[,1])

[1] 2233

## Option 2: Compute a frequency table of the values column 1
table(yose_poi_merced_mat[,1])

FALSE  TRUE 
  487  2233 

Subset the POIs that fall within the Upper Merced Subbasin

To extract the POIs in the Upper Merced subbasin, we can feed the first column of yose_poi_merced_mat into filter() (which expects TRUE/FALSE values):

## Extract the points that intersect the subbasin to a new sf object
merced_poi_utm <- yose_poi_utm %>% 
  filter(yose_poi_merced_mat[,1])

Plot the Intersection

Plot to visually verify the results:

## Plot
tm_shape(merced_huc8_utm_sf) +
  tm_polygons(col = "khaki") +
tm_shape(yose_poi_utm) +
  tm_dots(size = 0.1, col = "gray30") +
tm_shape(merced_poi_utm) +
  tm_dots(size = 0.1, col = "dodgerblue")

CHALLENGE: How many YNP POIs fall witin the Upper Tuolumne HUC-8 subbasin?

Answer

uptuol_huc8_utm_sf <- yose_huc8_sf %>% 
  filter(HUC_8_NAME == "UPPER_TUOLUMNE") %>% 
  st_transform(st_crs(yose_poi_utm))

## Plot to verify
tm_shape(uptuol_huc8_utm_sf) + 
  tm_polygons() +
tm_shape(yose_poi_utm) +
  tm_dots()

## Count the intersecting points
yose_poi_utm %>% 
  st_intersection(uptuol_huc8_utm_sf) %>% 
  st_drop_geometry() %>% 
  count()

attribute variables are assumed to be spatially constant throughout all geometries

End

Congratulations, you’ve completed another Notebook!

To view your Notebook at HTML, save it (again), then click the ‘Preview’ button in the RStudio toolbar.