Skip to content

Instantly share code, notes, and snippets.

@benmarwick

benmarwick/computing-sota-az-boundaries.R

Last active March 21, 2025 00:14
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save benmarwick/f7daa88453c2c5853e4e20ccc8e53b19 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save benmarwick/f7daa88453c2c5853e4e20ccc8e53b19 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. benmarwick revised this gist Mar 25, 2024. 1 changed file with 118 additions and 88 deletions.
    206 changes: 118 additions & 88 deletions computing-sota-az-boundaries.R
    View file
    Original file line number Diff line number Diff line change
    @@ -7,9 +7,10 @@ library(tidyverse)
    library(ggrepel)
    library(ggspatial)

    # get all summits via the GeoJSON download: https://sotl.as/summits/W7W/KG
    # 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")),
    @@ -19,18 +20,17 @@ gjsf %>%

    # 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 radius
    # challenge here is to make the buffer big enough
    # to include the AZ, but not too big to include two summits
    # get a buffer zone of a certain area
    buffer_side_length <- 1e6 # 100 is a 10x10
    gjsf_elev_buf <-
    bSquare(gjsf_elev, 1e6)
    bSquare(gjsf_elev, buffer_side_length)

    gjsf_elev_buf_sq <- vector("list", nrow(gjsf_elev_buf))
    for(i in 1:nrow(gjsf_elev_buf)){
    @@ -66,114 +66,139 @@ for(i in 1:nrow(gjsf_elev_buf_sq_df)){
    #----------------------------------
    for(i in 1:nrow(gjsf_elev)){

    which_summit <- i

    bbx = gjsf_elev_buf_sq_bbx[[which_summit]]
    this_summit_point <- gjsf_elev_buf[i, ]

    # as sf object
    bbx_st <- st_as_sf(bbx,
    coords = c("X", "Y"),
    crs = st_crs(gjsf_elev))
    print(paste0("Starting work on the AZ for ", this_summit_point$id,
    "..."))

    # as polygon
    bbx_st_poly <-
    st_convex_hull(st_union(bbx_st))
    # extract the bounding box for just this summit
    bbx_st_poly <- gjsf_elev_buf_sq_df[i, ]

    # quick look
    # quick look at this bbx in context
    ggplot() +
    geom_sf(data = bbx_st_poly, colour = "red") +
    geom_sf(data = gjsf_elev, size = 0.1)
    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 = gjsf_elev[which_summit,], size = 0.1)

    # subset points in this bbox
    where_at <- gjsf_elev_buf[which_summit,]

    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,
    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,
    z = 14, # this is the max possible
    prj = "WGS84")

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

    # find area that is in AZ in meters. If the SOTA summit elev is
    # less than the max elev in the raster, use that, otherwise
    # use the max elev in the raster. Because some SOTA summit elev
    # are higher than the elevation in the rasters
    where_at_az <-
    where_at %>%
    mutate(bbx_ras_max = bbx_ras_max,
    az_lower_contour = ifelse(elev_m <= bbx_ras_max,
    elev_m - 25,
    bbx_ras_max - 25))

    # elevation for Boise Ridge W7W/KG-114 might be wrong
    # max in raster is 941m, but SOTA data is 969m
    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 < where_at_az$az_lower_contour] <- NA

    # get extent as polygon
    pr <- st_as_sf(as.polygons(rast(rast) > -Inf))

    # if there are multiple, we only want the one that
    # contains the summit point when we have multipolys,
    # we just want the one with the summit in it

    multipoly <- pr
    summit_point <- gjsf_elev[which_summit,]
    df_union_cast <- st_cast(multipoly, "POLYGON")

    poly_with_summit <-
    apply(st_intersects(df_union_cast,
    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 = where_at,
    aes(x, y)) +
    geom_text_repel(data = gjsf_elev[which_summit,],
    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()

    library(geojsonsf)
    geo <- sf_geojson(poly_with_summit)

    file_name <- paste0("output/", str_replace_all(where_at$id, "/|-", "_"),
    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))

    # show in the console where
    print(file_name)

    }

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



  2. benmarwick revised this gist Mar 18, 2024. 1 changed file with 84 additions and 16 deletions.
    100 changes: 84 additions & 16 deletions computing-sota-az-boundaries.R
    View file
    Original file line number Diff line number Diff line change
    @@ -4,8 +4,10 @@ library(tidyterra)
    library(ggplot2)
    library(sf)
    library(tidyverse)
    library(ggrepel)
    library(ggspatial)

    # get all summits via the GeoJSON download from sotl.as
    # get all summits via the GeoJSON download: https://sotl.as/summits/W7W/KG
    gjsf = st_read("W7W_KG.geojson")

    gjsf_elev <-
    @@ -25,10 +27,11 @@ bSquare <- function(x, a) {
    }

    # get a buffer zone of a certain radius
    # challenge here is to make the buffer big enough
    # to include the AZ, but not too big to include two summits
    gjsf_elev_buf <-
    bSquare(gjsf_elev, 5e5)
    bSquare(gjsf_elev, 1e6)

    # really make it square
    gjsf_elev_buf_sq <- vector("list", nrow(gjsf_elev_buf))
    for(i in 1:nrow(gjsf_elev_buf)){
    gjsf_elev_buf_sq[[i]] <-
    @@ -44,10 +47,11 @@ bind_rows(gjsf_elev_buf_sq)
    ggplot() +
    geom_sf(data = gjsf_elev_buf_sq_df,
    colour = "red", fill = NA) +
    coord_sf()
    coord_sf() +
    annotation_scale(location = "bl",
    width_hint = 0.5)

    # get list of bounding boxes coords
    # as input for getting the DSM raster
    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]] <-
    @@ -77,7 +81,7 @@ bbx = gjsf_elev_buf_sq_bbx[[which_summit]]

    # quick look
    ggplot() +
    geom_sf(data = bbx_st_poly, colour = "red") +
    geom_sf(data = bbx_st_poly, colour = "red") +
    geom_sf(data = gjsf_elev, size = 0.1)

    ggplot() +
    @@ -87,32 +91,62 @@ ggplot() +
    # subset points in this bbox
    where_at <- gjsf_elev_buf[which_summit,]

    # download the raster tiles around the summit that have elevation data
    # Get the raster data from AWS
    # Elevation units are in meters.
    bbx_ras =
    elevatr::get_elev_raster(bbx_st,
    src = "aws",
    clip = "bbox",
    verbose = TRUE,
    z = 12,
    # 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,
    prj = "WGS84")

    # find area that is in AZ in meters
    # get max elev from raster
    bbx_ras_max <- maxValue(bbx_ras)

    # find area that is in AZ in meters. If the SOTA summit elev is
    # less than the max elev in the raster, use that, otherwise
    # use the max elev in the raster. Because some SOTA summit elev
    # are higher than the elevation in the rasters
    where_at_az <-
    where_at %>%
    mutate(az_lower_contour = elev_m - 25)
    mutate(bbx_ras_max = bbx_ras_max,
    az_lower_contour = ifelse(elev_m <= bbx_ras_max,
    elev_m - 25,
    bbx_ras_max - 25))

    # elevation for Boise Ridge W7W/KG-114 might be wrong
    # max in raster is 941m, but SOTA data is 969m

    rast <- bbx_ras
    rast[rast < where_at_az$az_lower_contour] <- NA

    # get extent as polygon
    pr <- st_as_sf(as.polygons(rast(rast) > -Inf))

    library(ggrepel)
    library(ggspatial)
    # if there are multiple, we only want the one that
    # contains the summit point when we have multipolys,
    # we just want the one with the summit in it

    multipoly <- pr
    summit_point <- gjsf_elev[which_summit,]
    df_union_cast <- st_cast(multipoly, "POLYGON")

    poly_with_summit <-
    apply(st_intersects(df_union_cast,
    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 = pr,
    geom_sf(data = poly_with_summit,
    colour ="red",
    fill = NA) +
    geom_point(data = where_at,
    @@ -128,20 +162,54 @@ ggplot() +
    theme_minimal()

    library(geojsonsf)
    geo <- sf_geojson(pr)
    geo <- sf_geojson(poly_with_summit)

    file_name <- paste0("output/", str_replace_all(where_at$id, "/|-", "_"),
    ".geojson")

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


    # show in the console where
    print(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 = 1) %>%
    addCircleMarkers(data = gjsf_elev,
    radius = 2,
    weight = 1)



  3. benmarwick created this gist Mar 17, 2024.
    149 changes: 149 additions & 0 deletions computing-sota-az-boundaries.R
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,149 @@
    library(raster)
    library(terra)
    library(tidyterra)
    library(ggplot2)
    library(sf)
    library(tidyverse)

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

    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

    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 radius
    gjsf_elev_buf <-
    bSquare(gjsf_elev, 5e5)

    # really make it square
    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()

    # get list of bounding boxes coords
    # as input for getting the DSM raster
    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)){

    which_summit <- i

    bbx = gjsf_elev_buf_sq_bbx[[which_summit]]

    # as sf object
    bbx_st <- st_as_sf(bbx,
    coords = c("X", "Y"),
    crs = st_crs(gjsf_elev))

    # as polygon
    bbx_st_poly <-
    st_convex_hull(st_union(bbx_st))

    # quick look
    ggplot() +
    geom_sf(data = bbx_st_poly, colour = "red") +
    geom_sf(data = gjsf_elev, size = 0.1)

    ggplot() +
    geom_sf(data = bbx_st_poly, colour = "red") +
    geom_sf(data = gjsf_elev[which_summit,], size = 0.1)

    # subset points in this bbox
    where_at <- gjsf_elev_buf[which_summit,]

    # download the raster tiles around the summit that have elevation data
    bbx_ras =
    elevatr::get_elev_raster(bbx_st,
    src = "aws",
    clip = "bbox",
    verbose = TRUE,
    z = 12,
    prj = "WGS84")

    # find area that is in AZ in meters
    where_at_az <-
    where_at %>%
    mutate(az_lower_contour = elev_m - 25)

    rast <- bbx_ras
    rast[rast < where_at_az$az_lower_contour] <- NA

    # get extent as polygon
    pr <- st_as_sf(as.polygons(rast(rast) > -Inf))

    library(ggrepel)
    library(ggspatial)

    ggplot() +
    geom_spatraster(data = rast(bbx_ras)) +
    geom_sf(data = pr,
    colour ="red",
    fill = NA) +
    geom_point(data = where_at,
    aes(x, y)) +
    geom_text_repel(data = gjsf_elev[which_summit,],
    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()

    library(geojsonsf)
    geo <- sf_geojson(pr)

    file_name <- paste0("output/", str_replace_all(where_at$id, "/|-", "_"),
    ".geojson")

    geojsonio::geojson_write(pr,
    file = here::here(file_name))


    print(file_name)

    }