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

## Your answer here

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

Answer

## Your answer here

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

Answer

## Your answer here

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

Answer

## Your answer here

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

## Your answer here

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

## Your answer here

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.