1. Zonal statistics with continuous data

1.1 Requirements

library(rgee)
library(sf)
library(viridis)   
library(tidyverse)
library(patchwork)
library(mapview)
library(ggspatial)
ee_Initialize()

── rgee 1.1.2 ────────────────────────── earthengine-api 0.1.292 ── 
 ✓ user: not_defined
 ✓ Initializing Google Earth Engine:  DONE!
 ✓ Earth Engine account: users/ambarja 
─────────────────────────────────────────────────────────────────── 

1.2 Vector layer reading of Lima

lima <- st_read(
  "https://github.com/healthinnovation/sdb-gpkg/raw/main/lima_distritos.gpkg",
  quiet = TRUE) %>% 
  dplyr::select(NOMBDIST,UBIGEO)

1.3 Transformation of sf object to a feature collection

lima_ee <-  lima %>% 
  sf_as_ee()

1.4 Processing data with rgee

# 1DU(Dobson units = 4,4615 x 10 -4 mol/m 2 molecules)
# source:https://sacs.aeronomie.be/info/dobson.php

so2 <- ee$ImageCollection$Dataset$COPERNICUS_S5P_NRTI_L3_SO2$
  select("SO2_column_number_density_15km")$
  filter(ee$Filter$calendarRange(2021,2021,"year"))$
  filter(ee$Filter$calendarRange(1,12,"month"))$
  median()$
  divide(0.00044615)$
  rename("so2")

# Extracting data
so2_to_lima <- ee_extract(
  x = so2,
  y = lima_ee,
  fun = ee$Reducer$mean(),
  sf = TRUE
)

glimpse(so2_to_lima)

Rows: 50
Columns: 4
$ NOMBDIST <chr> "LURIGANCHO", "JESUS MARIA", "LIMA", "LINCE", "MIRAFLORES", "M…
$ UBIGEO   <chr> "150118", "150113", "150101", "150116", "150122", "150120", "1…
$ so2      <dbl> 0.0267882, 0.0207625, 0.0148749, 0.0263857, 0.0321065, 0.00614…
$ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-76.7268 -1..., MULTIPOLYGON (((-…

1.5 Exploration Lima map

mapview(
  so2_to_lima,
  zcol="so2",
  layer.name = "SO2"
  )

1.6 Mapping with ggplot and ggspatial

p1 <- so2_to_lima %>% 
  ggplot() + 
  geom_sf(aes(fill = so2),color = NA) +
  scale_fill_viridis("SO2 [DU]",option = "F",direction = -1) + 
  theme_minimal(8) + 
  labs(
    title = "SO2 emission in Lima-2021"
    ) +
  annotation_north_arrow(
    height = unit(0.3,"cm"),
    width = unit(0.3,"cm"),
    location = "tr"
    ) + 
  annotation_scale(line_width = 0.05,height = unit(0.1,"cm"))

# Violin plot
v1 <-  so2_to_lima %>% 
  st_set_geometry(NULL) %>%
  mutate(name = "Lima Metropolitana") %>% 
  ggplot(aes(x = name ,y = so2)) + 
  geom_violin(width=1, size=0.5, alpha = 0.8,trim=FALSE) +
  geom_jitter(color = "darkgrey", alpha = 0.8, cex = 0.5) +
  theme_minimal() +
  labs(x = "",
      caption = "Source:Sentinel-5P Near Real-Time Sulphur Dioxide"
      )

p1 + v1

1.7 Vector layer reading of San Diego

sandiego <- st_read(
  "https://github.com/healthinnovation/sdb-gpkg/raw/main/SanDiego_districts.gpkg",
  quiet = TRUE) %>% 
  dplyr::select(OBJECTID,PO_NAME)

1.8 Transformation of sf object to a feature collection

sandiego_ee <- sandiego %>% 
  sf_as_ee()

1.9 Processing data with rgee

so2_to_sandiego <- ee_extract(
  x = so2,
  y = sandiego_ee,
  fun = ee$Reducer$mean(),
  sf = TRUE
)
glimpse(so2_to_sandiego)

Rows: 108
Columns: 4
$ OBJECTID <int> 28113, 28114, 28115, 28116, 28117, 28118, 28119, 28120, 28121,…
$ PO_NAME  <chr> "Alpine", "Bonita", "Boulevard", "Campo", "Chula Vista", "Chul…
$ so2      <dbl> 0.0157558, 0.0149694, 0.0139448, 0.0070619, 0.0198098, 0.02865…
$ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-116.6075 3..., MULTIPOLYGON (((-…

1.10 Exploration San Diego Map

mapview(
  so2_to_sandiego,
  zcol="so2",
  layer.name = "SO2")

1.11 Mapping with ggplot and ggspatial

p2 <- so2_to_sandiego %>% 
  ggplot() + 
  geom_sf(aes(fill = so2),color = NA) +
  scale_fill_viridis("SO2 [DU]",option = "F",direction = -1) + 
  theme_minimal(8) + 
  labs(
    title = "SO2 emission in San Diego-2021"
  ) +
  annotation_north_arrow(
    height = unit(0.3,"cm"),
    width = unit(0.3,"cm"),
    location = "tr"
  ) + 
  annotation_scale(line_width = 0.05,height = unit(0.1,"cm"))

# Violin plot
v2 <-  so2_to_sandiego %>% 
  st_set_geometry(NULL) %>%
  mutate(name = "San Diego") %>% 
  ggplot(aes(x = name ,y = so2)) + 
  geom_violin(width=1, size=0.5, alpha = 0.8,trim=FALSE) +
  geom_jitter(color = "darkgrey", alpha = 0.8, cex = 0.5) +
  theme_minimal() +
  labs(
    x = "",
    caption = "Source:Sentinel-5P Near Real-Time Sulphur Dioxide")+
  theme(plot.caption = element_text(hjust = 1, face = "italic"))

p2 + v2