Author

Carolyn Koehn

Please install and load these packages:

Code 
# install.packages('rgbif')

library(rgbif)
library(sf)
Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
Code 
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code 
library(spData)

Create new spatial attributes

Code 
nz.sf <- nz %>% 
  mutate(area = st_area(nz))
head(nz.sf$area, 3)
Units: [m^2]
[1] 12890576439  4911565037 24588819863
Code 
# generate random points
random_long_lat <- 
  data.frame(
    long = sample(runif(2000, min = 1090144, max = 2089533), replace = F),
    lat = sample(runif(2000, min = 4748537, max = 6191874), replace = F)
  )

random_long_lat_sf <- random_long_lat %>% 
  st_as_sf(coords = c("long", "lat"), crs = st_crs(nz))

random_nz <- random_long_lat_sf[nz.sf,]

# code from slides
nz.df <- nz %>% 
mutate(counts = lengths(st_intersects(., random_nz)),
       area = st_area(nz),
       density = counts/area)
head(st_drop_geometry(nz.df[,7:10]))
  counts              area              density
1     20 12890576439 [m^2] 1.551521e-09 [1/m^2]
2      4  4911565037 [m^2] 8.144044e-10 [1/m^2]
3     40 24588819863 [m^2] 1.626756e-09 [1/m^2]
4     21 12271015945 [m^2] 1.711350e-09 [1/m^2]
5      9  8364554416 [m^2] 1.075969e-09 [1/m^2]
6     18 14242517871 [m^2] 1.263821e-09 [1/m^2]
Code 
canterbury = nz %>% filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]
co = filter(nz, grepl("Canter|Otag", Name))
st_distance(nz_height[1:3, ], co)
Units: [m]
          [,1]     [,2]
[1,] 123537.16 15497.72
[2,]  94282.77     0.00
[3,]  93018.56     0.00
Code 
library(tigris)
To enable caching of data, set `options(tigris_use_cache = TRUE)`
in your R script or .Rprofile.
Code 
ua <- urban_areas(cb = FALSE, progress_bar = FALSE) %>% 
  filter(., UATYP10 == "U") %>% 
  filter(., str_detect(NAME10, "ID")) %>% 
  st_transform(., crs=2163)
Retrieving data for the year 2022
Warning in CPL_crs_from_input(x): GDAL Message 1: CRS EPSG:2163 is deprecated.
Its non-deprecated replacement EPSG:9311 will be used instead. To use the
original CRS, set the OSR_USE_NON_DEPRECATED configuration option to NO.
Code 
#get index of nearest ID city
nearest <-  st_nearest_feature(ua)
#estimate distance
(dist = st_distance(ua, ua[nearest,], by_element=TRUE))
Units: [m]
[1]  61373.575  61373.575   1647.128   1647.128 136917.546 136917.546
Code 
ua$dist_to_neighbor <- dist

Topological Subsetting

Code 
ctby_height <-  nz_height[canterbury, ]

ctby_height_diffpred <- nz_height[canterbury, ,op=st_touches]
Code 
canterbury_height3 = nz_height %>%
  filter(st_intersects(x = ., y = canterbury, sparse = FALSE))
Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
ℹ Please use one dimensional logical vectors instead.
Code 
# to fix warning
canterbury_height3 = nz_height %>%
  filter(as.vector(st_intersects(x = ., y = canterbury, sparse = FALSE)))

Spatial Joins

Code 
set.seed(2018)
(bb = st_bbox(world)) # the world's bounds
      xmin       ymin       xmax       ymax 
-180.00000  -89.90000  179.99999   83.64513
Code 
#>   xmin   ymin   xmax   ymax 
#> -180.0  -89.9  180.0   83.6

random_df = data.frame(
  x = runif(n = 10, min = bb[1], max = bb[3]),
  y = runif(n = 10, min = bb[2], max = bb[4])
)
random_points <- random_df %>% 
  st_as_sf(coords = c("x", "y")) %>% # set coordinates
  st_set_crs("EPSG:4326") # set geographic CRS

random_joined = st_join(random_points, world["name_long"])
Code 
any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))
[1] FALSE
Code 
z = st_join(cycle_hire, cycle_hire_osm, st_is_within_distance, dist = 20)
nrow(cycle_hire)
[1] 742
Code 
nc <- st_read(system.file("shape/nc.shp", package="sf")) %>%
  st_transform(crs=4326)
Reading layer `nc' from data source 
  `C:\Users\carolynkoehn\AppData\Local\R\cache\R\renv\cache\v5\R-4.3\x86_64-w64-mingw32\sf\1.0-16\ad57b543f7c3fca05213ba78ff63df9b\sf\shape\nc.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 100 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS:  NAD27
Code 
bb <- st_bbox(nc)
nc_points <- data.frame(
  x = runif(n = 2, min = bb[1], max = bb[3]),
  y = runif(n = 2, min = bb[2], max = bb[4])
) %>%
  st_as_sf(., coords = c("x", "y"), crs=4326)
tr_buff <- st_buffer(nc_points, units::set_units(100, "km")) %>%
  st_transform(crs=4326)

intersect_pct <- st_intersection(nc, tr_buff) %>% 
   mutate(intersect_area = st_area(.)) %>%   # create new column with shape area
   dplyr::select(NAME, intersect_area) %>%   # only select columns needed to merge
   st_drop_geometry()
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Code 
nc <- mutate(nc, county_area = st_area(nc))

# Merge by county name
nc <- left_join(nc, intersect_pct)
Joining with `by = join_by(NAME)`
Code 
# Calculate coverage
nc <- nc %>% 
   mutate(coverage = as.numeric(intersect_area/county_area)) %>%
  mutate(coverage = ifelse(test = is.na(coverage),
                           yes = 0,
                           no = coverage))

plot(nc["coverage"])

Practice

What is the relationship between observations of golden eagles and distance to a main road?

Code 
idaho <- tigris::states(progress_bar = FALSE) %>%
  filter(NAME == "Idaho")
Retrieving data for the year 2021
Code 
# get a selection of 1000 observations from rgbif
gold_eagles_us <- occ_search(scientificName = "Aquila chrysaetos", 
                    country = "US",
                    hasCoordinate = TRUE,
                    limit=1000)

# This step has to be separate, it can't be piped into the next bit of code
gold_eagles_us <- gold_eagles_us$data

# convert to spatial object
gold_eagle_dat_sf <- gold_eagles_us %>%
  filter(!is.na(decimalLatitude) & !is.na(decimalLongitude)) %>%
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)

# get roads
roads <- tigris::primary_secondary_roads("ID", progress_bar = FALSE) %>%
  st_transform(crs=4326)
Retrieving data for the year 2022
Code 
# check plot
plot(st_geometry(idaho))
plot(st_geometry(roads), col="red", add=TRUE)
plot(st_geometry(gold_eagle_dat_sf), add=TRUE)

Code 
# spatial subset to points in Idaho
idaho <- st_transform(idaho, crs=4326)

gold_eagle_dat_sf_id <- gold_eagle_dat_sf[idaho, ]

# check plot
plot(st_geometry(idaho))
plot(st_geometry(roads), col="red", add=TRUE)
plot(st_geometry(gold_eagle_dat_sf_id), add=TRUE)

Code 
# get nearest road for each point
nearest_road <- st_nearest_feature(gold_eagle_dat_sf_id, roads)

# get distance to nearest road for each point
gold_eagle_dat_sf_id <- gold_eagle_dat_sf_id %>%
  mutate(dist_to_road = st_distance(., roads[nearest_road,], by_element = TRUE))

# now we can investigate
# for example...
hist(gold_eagle_dat_sf_id$dist_to_road)