carbon works! 🌱

Author: cboettig

app/app.R

@@ -81,8 +81,8 @@ Smaller areas will be faster to compute!  Zoom in further to show richness with
     card(
       card_header("Biodiversity"),
       # chat_ui("chat", placeholder = "hummingbirds"),
-      actionLink("get_richness", "🐦 richness hexes"),
-      actionLink("get_mean_richness", "🐦 zonal richness")
+      actionLink("get_richness", "🐦 GBIF species richness"),
+      actionLink("get_carbon", "🌱 vulnerable carbon"),
     ),
     card(
       card_header("Resolution"),
@@ -93,7 +93,8 @@ Smaller areas will be faster to compute!  Zoom in further to show richness with
         max = MAXZOOM,
         value = 2,
         step = 1
-      )
+      ),
+      input_switch("show_hexes", "show hexes", value = FALSE)
     ),
     br(),
     input_switch("toggle_natgeo", "natgeo", value = TRUE),
@@ -128,7 +129,7 @@ server <- function(input, output, session) {
 
   taxa_selections <- taxonomicSelectorServer("taxa_selector")
 
-  # React to changes in taxonomic selections
+  # Update taxa selector tool to drill down on taxonomic selections
   observe({
     selections <- taxa_selections$selections()
     print(paste(
@@ -137,34 +138,9 @@ server <- function(input, output, session) {
     ))
     taxa_filter(selections)
   })
-  # Waterfall strategy to determine feature selection:
-  get_active_feature <- function(input) {
-    gdf <- active_feature()
-
-    if (is_empty(gdf)) {
-      print("No feature selected, checking for drawing")
-      gdf <- get_drawn_features(maplibre_proxy("map"))
-    }
-    if (is_empty(gdf)) {
-      print("No drawing found, checking geocoder")
-      gdf <- geocoder_to_gdf(input$map_geocoder)
-    }
-    if (is_empty(gdf)) {
-      print("No geocoder, getting current bbox")
-      bbox <- input$map_bbox
-      print(bbox)
-      if (!is.null(bbox)) {
-        gdf <- get_polygon_bbox(bbox)
-      }
-    }
-
-    if (is_empty(gdf)) {
-      warning("No selection found, using default!")
-      return(spData::us_states)
-    }
 
-    gdf
-  }
+  # Waterfall strategy to determine feature selection:
+  # FIXME select all (visible) protected areas? select protected areas by filtering them?
 
   # Set up the map:
   output$map <- renderMaplibre({
@@ -195,6 +171,7 @@ server <- function(input, output, session) {
     m <- m |> add_hillshade(visibility = "none")
 
     m <- m |> richness_layer() # default richness layer
+    m <- m |> carbon_layer() # default carbon layer
 
     m |> add_countries()
   })
@@ -239,30 +216,6 @@ server <- function(input, output, session) {
     proxy |> mapgl::set_filter(input$layer_selection, layer_filter())
   })
 
-  # Observe chat input
-  observeEvent(input$chat_user_input, {
-    taxa_selected <- txt_to_taxa(input$chat_user_input)
-
-    resp <- bot_response(taxa_selected, input$resolution)
-    print(resp)
-    chat_append("chat", resp)
-
-    # optionally - store the selection as global variable for future reactions
-    taxa_filter(taxa_selected)
-
-    # we can react right away, computing richness and updating map
-    gdf <- get_richness(
-      poly = get_active_feature(input),
-      zoom = as.integer(input$resolution),
-      taxa_selections = taxa_selected
-    )
-
-    chat_clear("chat")
-
-    maplibre_proxy("map") |>
-      set_source("richness", gdf)
-  })
-
   # Zoom into selected feature, move down a layer, show resulting child features
   observeEvent(input$map_feature_click, {
     x <- input$map_feature_click
@@ -316,35 +269,53 @@ server <- function(input, output, session) {
   })
 
   observeEvent(input$get_richness, {
-    gdf <- get_richness(
-      poly = get_active_feature(input),
-      zoom = as.integer(input$resolution),
-      taxa_selections = taxa_filter()
-    )
+    poly <- get_active_feature(active_feature(), input)
+
+    if (input$show_hexes) {
+      gdf <- get_richness(
+        poly = poly,
+        zoom = as.integer(input$resolution),
+        taxa_selections = taxa_filter()
+      )
+    } else {
+      layer <- layer_config[[input$layer_selection]]$parent_layer
+      child_poly <- child_polygons(poly, layer, layer_config)
+      gdf <- get_zonal_richness(
+        child_poly,
+        zoom = as.integer(input$resolution),
+        taxa_selections = taxa_filter()
+      )
+    }
 
     print(gdf)
     maplibre_proxy("map") |>
       set_source("richness", gdf)
   })
 
-  observeEvent(input$get_mean_richness, {
-    poly <- get_active_feature(input)
-    layer <- layer_config[[input$layer_selection]]$parent_layer
-    child_poly <- child_polygons(poly, layer, layer_config)
-    gdf <- get_zonal_richness(child_poly, zoom = as.integer(input$resolution))
+  observeEvent(input$get_carbon, {
+    poly <- get_active_feature(active_feature(), input)
+
+    if (input$show_hexes) {
+      gdf <- get_carbon(
+        poly = poly,
+        zoom = as.integer(input$resolution)
+      )
+    } else {
+      layer <- layer_config[[input$layer_selection]]$parent_layer
+      child_poly <- child_polygons(poly, layer, layer_config)
+      gdf <- get_mean_carbon(child_poly, zoom = as.integer(input$resolution))
+    }
 
-    print("mapping zonal richness...")
-    maplibre_proxy("map") |> set_source("richness", gdf)
+    print(gdf)
+    maplibre_proxy("map") |>
+      set_source("carbon", gdf)
   })
 
+  # Layer selection tools
   observeEvent(input$clear_filters, {
     maplibre_proxy("map") |> set_filter(input$layer_selection, NULL)
   })
 
-  observeEvent(input$clear_richness, {
-    maplibre_proxy("map") |> set_source("richness", "https://example.com")
-  })
-
   observeEvent(input$clear_area, {
     active_feature(NULL)
     layer_filter(NULL)
@@ -405,6 +376,30 @@ server <- function(input, output, session) {
         set_style(input$basemap)
     }
   })
+
+  # Observe chat input
+  observeEvent(input$chat_user_input, {
+    taxa_selected <- txt_to_taxa(input$chat_user_input)
+
+    resp <- bot_response(taxa_selected, input$resolution)
+    print(resp)
+    chat_append("chat", resp)
+
+    # optionally - store the selection as global variable for future reactions
+    taxa_filter(taxa_selected)
+
+    # we can react right away, computing richness and updating map
+    gdf <- get_richness(
+      poly = get_active_feature(active_feature(), input),
+      zoom = as.integer(input$resolution),
+      taxa_selections = taxa_selected
+    )
+
+    chat_clear("chat")
+
+    maplibre_proxy("map") |>
+      set_source("richness", gdf)
+  })
 }
 
 shinyApp(ui, server)

app/data-layers.R

@@ -193,6 +193,27 @@ richness_layer <- function(
     )
 }
 
+carbon_layer <- function(
+  map,
+  gdf = "https://minio.carlboettiger.info",
+  n_stops = 7
+) {
+  map |>
+    mapgl::add_source("carbon_source", gdf) |>
+    mapgl::add_fill_extrusion_layer(
+      id = "carbon",
+      source = "carbon_source",
+      tooltip = concat("carbon:", mapgl::get_column("carbon")),
+      fill_extrusion_color = mapgl::interpolate(
+        column = "value",
+        values = seq(0, 1, length.out = n_stops),
+        stops = viridisLite::viridis(n_stops, option = "viridis")
+      ),
+      fill_extrusion_height = list("*", 10000, list("get", "value")),
+      fill_extrusion_opacity = 0.7
+    )
+}
+
 
 add_richness_2d <- function(map, gdf) {
   map |>

app/tools-carbon.R

@@ -41,10 +41,10 @@ get_carbon <- function(
   poly,
   zoom = 8L,
   id_column = "id",
-  max_features = getOption("shiny_max_features", 20000L),
   warning = TRUE,
   verbose = TRUE,
-  server = Sys.getenv("AWS_S3_ENDPOINT", "minio.carlboettiger.info")
+  server = Sys.getenv("AWS_S3_ENDPOINT", "minio.carlboettiger.info"),
+  bucket = "public-data/cache/gbif-app"
 ) {
   duckdbfs::load_h3()
 
@@ -55,25 +55,46 @@ get_carbon <- function(
   ## This operation is maybe always fast enough not to cache?
   carbon <- open_carbon_region(poly_hexed, server) |>
     dplyr::group_by(h3id) |>
-    dplyr::summarise(carbon = mean(carbon))
-
-  # in-memory gdf will crash above a certain number of hexes
-  if (warning) {
-    n_features <- carbon |> count() |> pull(n)
-    print(paste("computed", n_features, "hexes"))
-    if (n_features > max_features) {
-      warning(paste("returning only first", max_features, "of", n_features))
-    }
-  }
+    dplyr::summarise(carbon = mean(carbon)) |>
+    dplyr::mutate(value = carbon / max(carbon)) # normalize for color scale
 
   carbon <- carbon |>
-    head(max_features) |> # max number of features
-    dplyr::mutate(geom = h3_cell_to_boundary_wkt(h3id))
+    dplyr::mutate(geom = ST_GeomFromText(h3_cell_to_boundary_wkt(h3id)))
 
-  hash <- digest::digest(list(gdf, zoom, id_column, label))
+  ## this part should be separate? Or be included in cache logic.
+  label <- "carbon"
+  hash <- digest::digest(list(carbon, zoom, id_column, label))
   s3 <- glue::glue("s3://{bucket}/{label}/{hash}.geojson")
-  duckdbfs::to_geojson(gdf, s3, as_http = TRUE)
+  duckdbfs::to_geojson(carbon, s3, as_http = TRUE)
 }
 
 
-get_mean_carbon <- function() {}
+get_mean_carbon <- function(
+  poly,
+  zoom = 8L,
+  id_column = "id",
+  warning = TRUE,
+  verbose = TRUE,
+  server = Sys.getenv("AWS_S3_ENDPOINT", "minio.carlboettiger.info"),
+  bucket = "public-data/cache/gbif-app"
+) {
+  # get_h3_aoi is self-caching, shared across metrics
+  poly_hexed_url <- get_h3_aoi(poly, precision = zoom, keep_cols = id_column)
+  poly_hexed <- duckdbfs::open_dataset(poly_hexed_url, recursive = FALSE)
+
+  ## This operation is maybe always fast enough not to cache?
+  carbon <- open_carbon_region(poly_hexed, server) |>
+    dplyr::group_by(.data[[id_column]]) |>
+    dplyr::summarise(carbon = mean(carbon)) |>
+    dplyr::mutate(value = carbon / max(carbon)) # normalize for color scale
+
+  gdf <- poly |>
+    dplyr::select(dplyr::all_of(id_column), geometry) |>
+    dplyr::inner_join(carbon, by = id_column) |>
+    rename(geom = "geometry")
+
+  label <- "carbon"
+  hash <- digest::digest(list(gdf, zoom, id_column, label))
+  s3 <- glue::glue("s3://{bucket}/{label}/{hash}.geojson")
+  duckdbfs::to_geojson(gdf, s3, as_http = TRUE)
+}