Spatial Queries - Sample Park Vegetation with Random Points

In this Notebook we’ll use geoprocessing functions from sf to:

• generate 500 randomly distributed sample points within YNP
• identify the vegetation class for sample point and save it in the attribute table
  
• find all sample points within 2km of each campground, and compute their distances to the campground

Setup

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

library(sf)
library(tibble)
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)

Import the Park Boundary

First, we import the boundary:

## Define a convenience variable for UTM Zone 11
epsg_utm11n_nad83 <- 26911

## Import the YNP border
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

Generate Random Points

Next, generate uniformly distributed random numbers between the x and y ranges of the park’s extent (i.e., bounding box). We’ll generate more points than we actually needed knowing that some of them will fall outside the park boundary.

n <- 500

## Get the Extent of the park boundary
yose_ext <- st_bbox(yose_bnd_utm)
yose_ext

     xmin      ymin      xmax      ymax 
 246310.4 4152875.3  306742.6 4229525.8 

## Generate uniformly distributed random numbers
xs <- runif(n * 3, min = yose_ext[1], max = yose_ext[3])
ys <- runif(n * 3, min = yose_ext[2], max = yose_ext[4])

Create a sf object from the random points

rand_pts_sf <- st_as_sf(data.frame(lon = xs, lat = ys, ptid = 1:length(xs)), 
                     coords = c("lon", "lat"), 
                     crs = st_crs(yose_bnd_utm))

## Plot to verify
tm_shape(rand_pts_sf) + 
  tm_dots(col = "dimgray") +
tm_shape(yose_bnd_utm) +
  tm_borders(col = "red", lwd = 2)

Next, test which points fall within the park boundary, and save those to a new object.

## Test which points are within the YNP boundary
in_yose_mat <- st_within(rand_pts_sf, yose_bnd_utm, sparse = FALSE)
head(in_yose_mat)

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

## Save those in the park to a new sf object
pts_in_park_sf <- rand_pts_sf %>% filter(in_yose_mat[,1])
nrow(pts_in_park_sf)

[1] 973

Plot the results:

tm_shape(yose_bnd_utm) + 
  tm_borders(col="palegreen4", lwd = 3) +
tm_shape(pts_in_park_sf) + 
  tm_dots(col="dimgray")

Lastly, all we need is 500 points, so take a random sample:

## Take a random sample of the points 
yose_sample_pts_utm <- pts_in_park_sf %>% 
  sample_n(n)

## Plot
tm_shape(yose_bnd_utm) + 
  tm_borders(col="palegreen4", lwd = 4) +
tm_shape(yose_sample_pts_utm) + 
  tm_dots(col="dimgray")

Find and Record the Vegetation Class for Each Sample Point

Next we’ll use a spatial join to do two things at once: 1) find the vegetation type for each sample point, and ii) save it in the attribute table.

First import the vegetation class layer and the legend:

yose_veg37_utm <- st_read(dsn="./data", layer="veg37") %>% 
  st_transform(epsg_utm11n_nad83) %>% 
  select(DOMINANT, ALLIANCE, RIM_CODE)

Reading layer `veg37' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data' using driver `ESRI Shapefile'
Simple feature collection with 6512 features and 26 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 246268.6 ymin: 4152706 xmax: 306269.9 ymax: 4229672
Projected CRS: NAD83 / UTM zone 11N

## Import the legend for the vegetation layer (saved separately)
yose_veg37_legend_df <- foreign::read.dbf("./data/veg37_alliances.dbf")

## Join the legend to the sf data frame on the 'ALLIANCE' column
yose_veg37leg_utm <- yose_veg37_utm %>% 
  left_join(yose_veg37_legend_df, by = "ALLIANCE")

yose_veg37leg_utm %>% slice(1:10) %>% as_tibble()

Plot the vegetation layer:

tmap_options(max.categories = 31)

tm_shape(yose_bnd_utm) + 
  tm_borders(col="palegreen4", lwd = 4) +
tm_shape(yose_veg37leg_utm) + 
  tm_fill("ALLIANCE") +
tm_layout(legend.outside = TRUE)

Now we’re ready to spatially join the sample points to the vegetation layer using st_join(). This will add the columns from the vegetation layer to the attribute table of the sample points.

yose_sample_pts_veg_utm <- yose_sample_pts_utm %>% 
  st_join(yose_veg37leg_utm)

yose_sample_pts_veg_utm %>% slice(1:10) %>% as_tibble()

CHALLENGE: Plot the sample points

Plot the sample points according to their vegetation class.

Answer

tm_shape(yose_bnd_utm) + 
  tm_borders(col="palegreen4", lwd = 4) +
tm_shape(yose_sample_pts_veg_utm) + 
  tm_dots("ALLIANCE", size = 0.3) +
tm_layout(legend.outside = TRUE)

Find All Sample Points within 2km of Each Campground

First we import the campgrounds:

## Import the campground
yose_campgrounds_utm <- st_read("./data", layer="yose_poi") %>% 
  st_transform(epsg_utm11n_nad83) %>% 
  filter(POITYPE == 'Campground') %>% 
  select(POINAME)

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

yose_campgrounds_utm

Simple feature collection with 15 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 247625.8 ymin: 4158885 xmax: 292815.1 ymax: 4204389
Projected CRS: NAD83 / UTM zone 11N
First 10 features:
                       POINAME                 geometry
1    Hodgdon Meadow Campground POINT (247625.8 4187296)
2        Crane Flat Campground POINT (253315.9 4181287)
3     Tamarack Flat Campground POINT (258936.4 4181679)
4        White Wolf Campground POINT (267142.6 4194768)
5    Yosemite Creek Campground POINT (271702.9 4189756)
6    Porcupine Flat Campground POINT (273987.5 4187789)
7  Tuolumne Meadows Campground POINT (292815.1 4194103)
8  Bridalveil Creek Campground POINT (268815.7 4171688)
9            Wawona Campground POINT (263419.9 4158885)
10           Camp 4 Campground   POINT (270612 4180317)

To identify the sample points within 2km of each campground, we could construct a buffer around each campground and take the intersection. But there’s an easier way - we can do the spatial query directly using st_is_within_distance().

The number of sample points near each campground will vary. Hence we’ll save the results in a list (sparse = TRUE):

samp_pts_near_each_campgrnd_lst <- yose_campgrounds_utm  %>% 
  st_is_within_distance(yose_sample_pts_utm, 2000, sparse = TRUE)

samp_pts_near_each_campgrnd_lst

Sparse geometry binary predicate list of length 15, where the predicate was `is_within_distance'
first 10 elements:
 1: (empty)
 2: 362
 3: 70
 4: 333
 5: (empty)
 6: 211, 375, 383, 416
 7: 97, 125
 8: 39, 102, 111, 143, 393
 9: 14
 10: 18, 55, 75, 327, 346, 358, 445

Compute Distances Between each Campground and its Nearby Sample Points

Finally, We’ll compute the distance between each campground and the sample points that lie within 2km.

To do this, we’ll loop through the campgrounds, and for each one we’ll use st_distance() to find the distance between the campground and the nearby points (whose row numbers are saved in samp_pts_near_each_campgrnd_lst) . To complile all the results, we’ll append the distances to a data frame each iteration of the loop.

## Initialize campgrnd_samppts_dist_df to NULL
## (we'll use it below to compile and save the results of a loop)
campgrnd_samppts_dist_df <- NULL

for (i in 1:nrow(yose_campgrounds_utm)) {
  
  if (length(samp_pts_near_each_campgrnd_lst[[i]]) > 0) {
    dist_mat <- st_distance(x = yose_campgrounds_utm %>% slice(i),
                          y = yose_sample_pts_utm %>% 
                            slice(samp_pts_near_each_campgrnd_lst[[i]]))
    
    ## Create a data frame with the campground rown number, the sample point row numbers, and the distances
    dist2samppts_df <- data.frame(campgrnd_idx = i,
                                  samp_pt_idx = samp_pts_near_each_campgrnd_lst[[i]],
                                  dist = as.numeric(dist_mat[1,]))
    
    ## Add these rows to campgrnd_samppts_dist_df
    campgrnd_samppts_dist_df <- campgrnd_samppts_dist_df %>% 
      bind_rows(dist2samppts_df)
    
  }

}

## View results
campgrnd_samppts_dist_df

End

Congratulations, you’ve completed another Notebook!

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