Compute polygons for SOTA activation zone boundaries for the W7W Association

computing-sota-az-boundaries.R

library(raster)
library(terra)
library(tidyterra)
library(ggplot2)
library(sf)
library(tidyverse)
library(ggrepel)
library(ggspatial)

# get all summits via the GeoJSON download
gjsf = st_read("W7W_KG.geojson")

# make sf data frame and get coords into cols we can use
gjsf_elev <- 
gjsf %>% 
  mutate(elev_m = parse_number(str_extract(title, ",\\s*(.*?)\\s*m")),
         elev_ft = elev_m * 3.2808399) %>% 
  mutate(x = st_coordinates(geometry)[,"X"],
         y = st_coordinates(geometry)[,"Y"])

# make a single square buffer for each point
# Function to create the square buffers
# https://stackoverflow.com/a/70372149/1036500
bSquare <- function(x, a) {
  a <- sqrt(a)/2
  return( sf::st_buffer(x, dist = a, nQuadSegs=1, 
                        endCapStyle = "SQUARE") )
}

# get a buffer zone of a certain area
buffer_side_length <- 1e6 # 100 is a 10x10
gjsf_elev_buf <- 
bSquare(gjsf_elev, buffer_side_length)

gjsf_elev_buf_sq <- vector("list", nrow(gjsf_elev_buf))
for(i in 1:nrow(gjsf_elev_buf)){
  gjsf_elev_buf_sq[[i]] <- 
    st_bbox(gjsf_elev_buf[i,]) %>%
    st_as_sfc() %>%
    st_as_sf()
}

gjsf_elev_buf_sq_df <- 
bind_rows(gjsf_elev_buf_sq)

# take a look at the buffer zones
ggplot() + 
  geom_sf(data = gjsf_elev_buf_sq_df,
          colour = "red", fill = NA) +
  coord_sf() +
  annotation_scale(location = "bl", 
                   width_hint = 0.5) 

# get list of bounding boxes coords
gjsf_elev_buf_sq_bbx <- vector("list", nrow(gjsf_elev_buf_sq_df))
for(i in 1:nrow(gjsf_elev_buf_sq_df)){
  gjsf_elev_buf_sq_bbx[[i]] <- 
    gjsf_elev_buf_sq_df$x[i] %>% 
    st_coordinates() %>% 
    as.data.frame() %>% 
    select(X, Y)
}

# loop over each bounding box, get the DSM raster 
# and find the AZ polygon and save as geoJSON
#----------------------------------
for(i in 1:nrow(gjsf_elev)){
  
this_summit_point <- gjsf_elev_buf[i, ]
  
print(paste0("Starting work on the AZ for ", this_summit_point$id,
               "..."))
  
# extract the bounding box for just this summit
bbx_st_poly <- gjsf_elev_buf_sq_df[i, ]
  
  # quick look at this bbx in context
ggplot()   +  
    geom_sf(data = bbx_st_poly, colour = "red") +
    geom_sf(data = gjsf_elev, size = 0.1) +
  theme_minimal() +
  coord_sf() +
  annotation_scale(location = "bl", 
                   width_hint = 0.5) 

# quick look at this bbx by itself
ggplot()   +  
  geom_sf(data = bbx_st_poly, colour = "red") +
  geom_sf(data = this_summit_point, size = 0.1) +
  theme_minimal() +
  coord_sf() +
  annotation_scale(location = "bl", 
                   width_hint = 0.5) 

#  Get the raster data from AWS
# Elevation units are in meters.
bbx_ras =
  elevatr::get_elev_raster(bbx_st_poly %>% 
                             st_cast("POINT"), 
                           src = "aws",
                           clip = "bbox", 
                           verbose = TRUE, 
                           # resolution info is here 
                           # https://github.com/tilezen/joerd/blob/master/docs/data-sources.md#what-is-the-ground-resolution
                           # around 5 m / pixel at z = 14 ?
                           z = 14, # this is the max possible
                           prj =  "WGS84")

# get max elev from raster
bbx_ras_max <- maxValue(bbx_ras)

print(paste0("Max elevation in raster is ", bbx_ras_max,
             " and summit elev is ", this_summit_point$elev_m))

# find area that is in AZ in meters. If the SOTA summit elev is 
# greater than the max elev in the raster, this is an error, so
# let's handle it. e.g. elevation for Boise Ridge W7W/KG-114 
# might be wrong: max in raster is 941m, but SOTA data is 969m

# define elevation contour that bounds the AZ
this_summit_point_az <- 
    this_summit_point %>% 
    mutate(bbx_ras_max = bbx_ras_max, # this is helpful for debugging
           az_lower_contour = ifelse(elev_m <= bbx_ras_max,
                                     elev_m - 25,         # AZ is area -25m elevation from summit
                                     bbx_ras_max - 25 ))   # SOTA summit data does not always match raster data

# subset the summit point that is in this bbox
rast <- bbx_ras
  rast[rast < this_summit_point_az$az_lower_contour] <- NA
  
# get extent of the AZ raster as polygon
  az_poly <- st_as_sf(as.polygons(rast(rast) > -Inf))
  
  # if there are multiple polygons, we only want the one that 
  # contains the summit point when we have multipolys, 
  # we just want the one with the summit in it
  
  df_union_cast <- st_cast(az_poly, "POLYGON")
  
  poly_with_summit <- 
    apply(st_intersects(df_union_cast, 
                        this_summit_point, 
                        sparse = FALSE), 2, 
          function(col) { 
            df_union_cast[which(col), ]
          })[[1]]
  
  # take a look at the summit and AZ
  ggplot() +
    geom_spatraster(data = rast(bbx_ras)) +
    geom_sf(data = poly_with_summit,
            colour ="red",
            fill = NA) +
    geom_point(data = this_summit_point,
               aes(x, y)) + 
    geom_text_repel(data = this_summit_point,
                    aes( x, y,
                         label = name),
                    bg.color = "white",
                    bg.r = 0.1) +
    scale_fill_viridis_c(na.value = "white") +
    annotation_scale(location = "bl", width_hint = 0.5) +
    coord_sf() +
    theme_minimal()

  # what is the longest linear dimension of this AZ polygon?
  # if it's the same as the bbox dimensions, then our bbox is 
  # probably too small and we need to get a bigger one
  poly_with_summit_max_linear_dim <- 
  poly_with_summit %>% 
        st_cast('MULTIPOINT') %>% 
        st_cast('POINT') %>% 
        st_distance %>% 
        max()
  
  class(poly_with_summit_max_linear_dim) <- "numeric"
  
  if(poly_with_summit_max_linear_dim < sqrt(buffer_side_length)) {
    
    print(paste0("OK: the max linear dimenstion of the AZ computed for ",
                 this_summit_point$id,
                 " is smaller than our bounding box"))
    
    print(paste0("The AZ for ", this_summit_point$id,
                 " was successfully computed"))
  } else {
     
    print(paste0("Not OK: the max linear dimenstion of the AZ computed for ",
                 this_summit_point$id,
                 " is bigger than our bounding box. Try a bigger bounding box."))
    
  }
  
# write AZ polygon to a GeoJSON file
file_name <- paste0("output/", str_replace_all(this_summit_point$id, "/|-", "_"), 
                    ".geojson")

# export AZ polygon as a GeoJSON file
geojsonio::geojson_write(poly_with_summit, 
                         file = here::here(file_name))
}

#------------------------------------------------------------
# browse the results to see how it looks

library(tidyverse)
library(sf)

# import the GeoJSON files with the AZ polygons
gjson_files <- 
  list.files(path = "output",
             full.names = TRUE,
             recursive = TRUE)

az_files <- 
  map(gjson_files,
      st_read) 

names(az_files) <- str_remove_all(basename(gjson_files), 
                                  "output|.geojson")

summit_geometry <- 
  bind_rows(az_files, .id = "summit") %>% 
  select(summit, geometry) 

# plot an interactive map with a nice base layer
library(leaflet)

leaflet() %>% 
  addProviderTiles("OpenStreetMap") %>% 
  addPolygons(data = summit_geometry,
              color = "red", 
              weight = 5) %>% 
  addCircleMarkers(data = gjsf_elev %>% 
                     mutate(id = str_replace_all(id, "/|-", "_")) %>% 
                     filter(id %in% summit_geometry$summit ),
                   radius = 2,
                   weight = 1,
                   label = ~id,
                   labelOptions = labelOptions(noHide = T, direction = "top",
                                               offset = c(0, -15)))