Module12

lucyengar@gmail.com

Module 11: Spatial Data

Using R for spatial data and plots.

Part 1

library(sf)

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(ggplot2)
library(viridis)    

## Loading required package: viridisLite

library(RColorBrewer)


countries <- st_read("countries.shp", quiet = TRUE)

if (is.na(st_crs(countries)))          
  st_crs(countries) <- 4326

make_choro <- function(data, var, trans = identity, palette = "viridis",
                       title = NULL, legend = NULL) {
  ggplot(data) +
    geom_sf(aes(fill = trans(.data[[var]])), colour = NA) +
    coord_sf(crs = st_crs(4326)) +
    scale_fill_viridis_c(option = palette, na.value = "grey80",
                         name = legend) +
    theme_minimal(base_size = 12) +
    labs(title = title)
}


gdp_map <- make_choro(countries,
                      var   = "gdp_md_est",
                      palette = "magma",
                      trans = log10,
                      title = "Median GDP (log10 scale)",
                      legend = "log10(GDP, million USD)")

pop_map <- make_choro(countries,
                      var   = "pop_est",
                      trans = sqrt,
                      title = "Population",
                      legend = "sqrt(Population)")
"Exercise 1"

## [1] "Exercise 1"

gdp_map

## Warning in FUN(X[[i]], ...): NaNs produced

pop_map

## Warning in trans(.data[["pop_est"]]): NaNs produced

Part 2

"Exercise 2"

## [1] "Exercise 2"

library(sf)
library(terra)

## terra 1.8.54

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:terra':
## 
##     intersect, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

b4 <- rast("LC08_044034_20170614_B4.tif")
b5 <- rast("LC08_044034_20170614_B5.tif")

ndvi <- (b5 - b4) / (b5 + b4)

places <- st_read("ca_places.shp", quiet = TRUE) |>
          st_transform(crs(ndvi))         

targets <- places |> 
           filter(NAME %in% c("Bethel Island", "Oakley"))

med <- extract(ndvi, vect(targets), fun = median, na.rm = TRUE)

data.frame(place = targets$NAME,
           median_NDVI = med[, 2])

##           place median_NDVI
## 1        Oakley   0.3005387
## 2 Bethel Island   0.4180438