Code for figure

  • Optional: examine the code and identify which parts of the code make each component of the figure.

Spatial continued

Data

districts <- st_read(system.file("extdata/districts.geojson", 
                                 package = "geospaar"))
farmers_sf <- read_csv(system.file("extdata/farmer_spatial.csv", 
                                   package = "geospaar"))
farmers_s <- st_as_sf(farmers, coords = c("x", "y")) %>% 
  st_set_crs(st_crs(districts))

roads <- st_read(system.file("extdata/roads.geojson", package = "geospaar")) %>% 
  st_transform(x = roads, crs = st_crs(districts))

pt <- st_point(x = c(28, -14))
pts <- st_multipoint(x = cbind(x = c(27.5, 28, 28.5), y = c(-14.5, -15, -15.5)))
sline <- st_linestring(cbind(x = c(27, 27.5, 28), y = c(-15, -15.5, -16)))
pol <- st_polygon(list(cbind(x = c(26.5, 27.5, 27, 26, 26.5), 
                             y = c(-15.5, -16.5, -17, -16, -15.5))))

Plotting (ggplot)

  • Use geom_sf()
ggplot(districts) + geom_sf() +
  geom_sf(data = roads, col = "red") + 
  geom_sf(data = farmers_sf, col = "blue")

Area and Length

  • st_area(), st_length. Default units are in meters
  • How many square m in a square km?
dist_areas <- districts %>% st_area()
dist_areas
dist_areas %>% units::set_units("km^2")

Add area as column

districts <- districts %>% 
  mutate(area = st_area(.) %>% units::set_units("km^2") %>% as.numeric()) 
districts

Plot districts (exercise)

  • Use geom_sf() to plot area
  • Fill districts based on area column.
  • Use scale_fill_continuous , with type = 'viridis'
ggplot(districts) + geom_sf(aes(fill = area)) + 
  scale_fill_continuous(type = "viridis")

Compare areas

districts <- districts %>% 
  st_transform("ESRI:102022") %>% 
  mutate(area2 = as.numeric(st_area(.) / 10000)) %>% 
  st_transform(4326)

ggplot(districts) + 
  geom_histogram(aes(x = area - area2))
districts <- districts %>%
  mutate(areadiff = (area - area2) / area * 100)

ggplot(districts) +
  geom_point(aes(x = area2, y = areadiff))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>%
  ggplot() + geom_point(aes(x = area2, y = areadiff))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>% 
  summarize(mu = mean(areadiff))
districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>%
  ggplot() + geom_sf(aes(fill = areadiff))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>% 
  summarize(mu = mean(areadiff)) %>% 
  ggplot() + geom_sf(aes(fill = mu))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>% 
  st_drop_geometry() %>% 
  summarize(mu = mean(areadiff)) #%>% 
  # ggplot() + geom_sf(aes(fill = mu))

districts %>%
  mutate(areadiff = ((area - area2) / area) * 100) %>% 
  as_tibble() %>% 
  summarize(mu = mean(areadiff)) #%>%

Calculate length (exercise)

  • Use st_length() to add column ‘road_length’ to roads
  • ‘road_length’ should be in km
  • Use ggplot to plot districts with black outline (color), and empty fill (fill = NA)
  • Then plot roads based on length (use scale_color_continuous() )
  • Plot roads with line width 1 (lwd = 1)
  • Use theme_bw()
roads <- st_read(
  system.file("extdata/roads.geojson", package = "geospaar")
) 
roads <- roads %>% 
  mutate(road_length = as.numeric(units::set_units(st_length(.), "km"))) 
ggplot() +
  geom_sf(data = districts, color = "black", fill = NA) +
  geom_sf(data = roads, aes(color = road_length), lwd = 1) + 
  scale_color_viridis_c() + 
  theme_minimal()

ggplot steps

  • Start with ggplot(). You can include data if you’re only using one data set.
  • Add plot type (geom_point, geom_line, geom_sf, geom_histogram)
  • Add aesthetics in aes(). These are columns used for plotting.
  • For geom_point, geom_line you need ‘x’ and ‘y’ aesthetics.
  • For geom_sf, spatial context used for ‘x’ and ‘y’. But you might want ‘fill’ or ‘color’