Geoprocessing 1 - Geometric Measurements

In this Notebook we will use functions from sf that measure geometric properties (area, length, distance).

• compute the total area of YNP in meters² and miles²
  
• compute the area of watersheds
  
• compute the total length of roads
  
• compute the distance between campgrounds and cell towers

Setup

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

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

tmap mode set to plotting

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

## 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

Compute the area of the park

Find area with st_area():

yose_area <- yose_bnd_utm %>% st_area()
yose_area

3019918537 [m^2]

Find the area in miles².

## View in square miles
set_units(yose_area, mi^2)

1165.997 [mi^2]

Import the Watersheds

Next import the watershed boundaries.

yose_watersheds_utm <- st_read("./data/yose_watersheds.gpkg", layer="calw221") %>% 
  st_transform(epsg_utm11n_nad83) %>% 
  select(CALWNUM, HU, HRNAME, RBNAME, HUNAME, HANAME, CDFSPWNAME, CDFPWSNAME, HUC_8, HUC_8_NAME)

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 38 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 1383.82 ymin: -61442.93 xmax: 81596.71 ymax: 26405.66
Projected CRS: unnamed

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

tm_shape(yose_watersheds_utm) +
  tm_borders() +
tm_shape(yose_bnd_utm) +
  tm_borders(col="red", lwd=3)

CHALLENGE: Compute the watershed areas

Compute the watershed areas, and add them as a new column in the attribute table. Answer

## Compute the area in map units
yose_watersheds_areas_m2 <- st_area(yose_watersheds_utm)

## Add the area to the attribute table
yose_watersheds_utm <- yose_watersheds_utm %>% 
  mutate(area_m2 = yose_watersheds_areas_m2)

## Preview results
yose_watersheds_utm %>% slice(1:10) %>% as_tibble()

CHALLENGE: Find the largest watershed

Find the largest watershed, and report the area in acres. Answer

## Find the biggest watershed
yose_watersheds_utm %>% 
  st_drop_geometry() %>%     ## get rid of the geometry column - don't need it
  top_n(1, area_m2) 

## Compute the area in acres
yose_watersheds_utm %>% 
  st_drop_geometry() %>%     ## get rid of the geometry column - don't need it
  top_n(1, area_m2) %>% 
  pull(area_m2) %>% 
  set_units("acres")

20017.12 [acres]

Import the Roads

yose_roads_utm <- st_read("./data/yose_roads.gdb", "Yosemite_Roads") %>% 
  st_transform(epsg_utm11n_nad83) %>% 
  select(RDNAME, RTENUMBER, YOSE_Surface, YOSE_FIRE_ROAD, YOSE_INPARK, YOSE_Type)

Reading layer `Yosemite_Roads' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data\yose_roads.gdb' using driver `OpenFileGDB'
Simple feature collection with 823 features and 40 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 234658.1 ymin: 4139484 xmax: 324852.6 ymax: 4250252
Projected CRS: NAD83 / UTM zone 11N

Compute the length of roads and save it in the attribute table:

yose_roads_utm <- yose_roads_utm %>% 
  mutate(rd_length = st_length(yose_roads_utm))
yose_roads_utm %>% slice(1:10) %>% as_tibble()

CHALLENGE: Compute Road Lengths

Find the total length for each type of Road (see YOSE_Type) in miles. Answer

yose_roads_utm %>% 
  st_drop_geometry() %>% 
  group_by(YOSE_Type) %>% 
  summarise(length_miles = sum(rd_length) %>% set_units("mi")) %>% 
  arrange(desc(length_miles))

Import the Campgrounds and Cell Towers

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_celltwrs_utm <- st_read("./data/yose_communications.gdb", "Cell_Towers") %>% 
  st_transform(epsg_utm11n_nad83) 

Reading layer `Cell_Towers' from data source `D:\Workshops\R-Spatial\rspatial_mod\outputs\rspatial_data\data\yose_communications.gdb' using driver `OpenFileGDB'
Simple feature collection with 5 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 251532.4 ymin: 4158756 xmax: 293307.2 ymax: 4194328
Projected CRS: NAD83 / UTM zone 11N

Plot them:

tm_shape(yose_bnd_utm) +
  tm_borders() +
tm_shape(yose_campgrounds_utm) +
  tm_symbols(size = 0.3, col = "blue") +
tm_shape(yose_celltwrs_utm) +
  tm_symbols(size = 0.3, col = "red") 

Compute Distances Between Features - 1 sf data frame

st_distance() can be used to compute distances between features. If you pass it one sf object, it will return a symmetric distance matrix between all pairs of features.

Compute the distance between all pairs of campgrounds:

campgrounds_dist_mat <- st_distance(yose_campgrounds_utm)
campgrounds_dist_mat

Units: [m]
           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]     [,7]      [,8]     [,9]     [,10]      [,11]
 [1,]     0.000  8275.796 12628.738 20897.961 24202.415 26366.307 45699.01 26318.172 32505.98 24022.360 27053.5366
 [2,]  8275.796     0.000  5634.176 19310.626 20243.499 21669.864 41526.22 18231.678 24574.95 17323.278 20522.7859
 [3,] 12628.738  5634.176     0.000 15448.311 15106.794 16243.769 36084.79 14050.939 23230.46 11754.726 14931.6215
 [4,] 20897.961 19310.626 15448.311     0.000  6776.104  9775.471 25681.09 23140.503 36074.82 14860.874 15706.6234
 [5,] 24202.415 20243.499 15106.794  6776.104     0.000  3014.800 21555.06 18297.281 31962.27  9501.164  9443.6883
 [6,] 26366.307 21669.864 16243.769  9775.471  3014.800     0.000 19858.20 16911.114 30774.47  8198.298  7236.6861
 [7,] 45699.014 41526.218 36084.793 25681.086 21555.058 19858.200     0.00 32839.145 45873.18 26134.546 23327.6854
 [8,] 26318.172 18231.678 14050.939 23140.503 18297.281 16911.114 32839.14     0.000 13892.94  8814.694 10183.5588
 [9,] 32505.983 24574.955 23230.456 36074.816 31962.268 30774.468 45873.18 13892.936     0.00 22606.591 24037.3432
[10,] 24022.360 17323.278 11754.726 14860.874  9501.164  8198.298 26134.55  8814.694 22606.59     0.000  3222.2366
[11,] 27053.537 20522.786 14931.621 15706.623  9443.688  7236.686 23327.69 10183.559 24037.34  3222.237     0.0000
[12,] 27174.786 20589.170 15009.085 15972.557  9717.027  7503.431 23439.68  9979.152 23820.38  3268.848   273.4135
[13,] 27088.787 20464.194 14892.172 16067.174  9849.639  7671.851 23645.40  9768.727 23613.03  3140.976   440.3027
[14,] 27616.957 20889.582 15345.375 16750.368 10495.364  8244.331 23684.28  9487.701 23272.44  3604.361  1053.9344
[15,] 22073.083 24539.640 22864.646 11107.445 17786.365 20717.442 32873.41 33489.566 45540.25 25705.918 26791.2639
           [,12]      [,13]      [,14]    [,15]
 [1,] 27174.7857 27088.7873 27616.9569 22073.08
 [2,] 20589.1697 20464.1936 20889.5823 24539.64
 [3,] 15009.0854 14892.1717 15345.3754 22864.65
 [4,] 15972.5573 16067.1742 16750.3676 11107.45
 [5,]  9717.0274  9849.6393 10495.3635 17786.37
 [6,]  7503.4311  7671.8512  8244.3308 20717.44
 [7,] 23439.6793 23645.4050 23684.2839 32873.41
 [8,]  9979.1521  9768.7266  9487.7011 33489.57
 [9,] 23820.3830 23613.0283 23272.4414 45540.25
[10,]  3268.8477  3140.9757  3604.3606 25705.92
[11,]   273.4135   440.3027  1053.9344 26791.26
[12,]     0.0000   212.8079   782.1046 27054.79
[13,]   212.8079     0.0000   692.4943 27142.02
[14,]   782.1046   692.4943     0.0000 27829.09
[15,] 27054.7878 27142.0250 27829.0938     0.00

Compute Distances Between Features - 2 sf data frames

If you pass two sf objects to st_distance(), it will return a n x m matrix containing the distance between all pairs of features.

Compute the distance between all pairs of campgrounds and cell towers:

cgct_dist_mat <- st_distance(yose_campgrounds_utm, yose_celltwrs_utm)
cgct_dist_mat

Units: [m]
            [,1]      [,2]      [,3]      [,4]      [,5]
 [1,] 46219.4109 25884.856 26225.834 33538.244  5903.577
 [2,] 42063.8195 19210.673 19163.004 25491.865  2384.789
 [3,] 36624.3349 13645.182 13737.613 23773.991  7499.216
 [4,] 26168.3039 15613.133 17308.435 36061.840 19627.095
 [5,] 22082.7907  9689.048 11543.363 31666.279 21313.320
 [6,] 20396.4206  7834.371  9723.708 30321.762 22987.367
 [7,]   541.1383 24682.675 25995.385 44830.662 42783.486
 [8,] 33352.8435  9178.442  7376.002 13416.944 20585.520
 [9,] 46361.8785 23070.914 21252.438  1824.674 26769.154
[10,] 26671.4314  1890.457  2646.959 22220.557 19249.627
[11,] 23859.2190  1411.166  2821.810 23427.597 22413.196
[12,] 23970.1116  1403.475  2639.230 23199.118 22497.225
[13,] 24175.5888  1257.224  2426.427 22995.829 22383.935
[14,] 24211.1160  1733.602  2314.224 22611.632 22843.581
[15,] 33272.4448 26622.758 28245.217 45778.512 23754.044

CHALLENGE: Find the closest cell tower

For the first campground (row), which is the closest cell tower? Answer

Hint: use which.min()

which.min(cgct_dist_mat[1,])

[1] 5

## Or for call campgrounds:
apply(cgct_dist_mat, 1, which.min)

 [1] 5 5 5 2 2 2 1 3 4 2 2 2 2 2 5

End

Congratulations, you’ve completed the Notebook!

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