render both

boettiger-lab · Dec 5, 2023 · d1ef998 · d1ef998
1 parent f19dd4b
commit d1ef998
Show file tree

Hide file tree

Showing 10 changed files with 78 additions and 25 deletions.
diff --git a/_freeze/contents/ace/execute-results/html.json b/_freeze/contents/ace/execute-results/html.json
@@ -1,7 +1,7 @@
 {
-  "hash": "77107a53791af3118af7d6dc5a464fab",
+  "hash": "997de147c6ee92af63b756f4c50126d6",
   "result": {
-    "markdown": "---\n\n---\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(stars)\nlibrary(rstac)\nlibrary(gdalcubes)\nlibrary(tmap)\nlibrary(dplyr)\ntmap_mode(\"plot\")\n\ngdalcubes::gdalcubes_options(parallel=24*2)\n```\n:::\n\n\nCDFW makes the [ACE GIS Data](https://wildlife.ca.gov/Data/Analysis/Ace) freely available, but in one large zip archive.\n\nHere I'll use a public mirror on a Berkeley-based server for a more cloud-native access pattern so we don't have to download the whole thing.  We also index this in a simple [stac catalog](https://radiantearth.github.io/stac-browser/#/external/raw.githubusercontent.com/schmidtDSE/biodiversity-catalog/main/stac/v1/collections/cdfw_ace/summary.json?.asset=asset-ds2715)\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nurl <- \"/vsicurl/https://minio.carlboettiger.info/public-biodiversity/ACE_summary/ds2715.gdb\"\n```\n:::\n\n\nCDFW uses hexes!  Note that these pre-date the modern [H3 hex tiling system](https://h3geo.org/), so do not provide the magic of [hierarchical (zoom level) tiles](https://h3geo.org/docs/highlights/indexing)\n\n\n::: {.cell}\n\n```{.r .cell-code}\nhexes <- st_read(url) \n```\n\n::: {.cell-output .cell-output-stdout}\n```\nReading layer `ds2715' from data source \n  `/vsicurl/https://minio.carlboettiger.info/public-biodiversity/ACE_summary/ds2715.gdb' \n  using driver `OpenFileGDB'\nSimple feature collection with 63890 features and 17 fields\nGeometry type: MULTIPOLYGON\nDimension:     XY\nBounding box:  xmin: -373987.9 ymin: -604495.8 xmax: 540034.3 ymax: 450023.2\nProjected CRS: NAD83 / California Albers\n```\n:::\n\n```{.r .cell-code}\ntm_shape(hexes) + tm_polygons(fill = \"RwiRankEco\", col=NA)\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-4-1.png){width=672}\n:::\n:::\n\n\nFind Biodiversity Intactness Index COG tiles from Planetary Computer using STAC search:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nbox <- st_bbox(st_transform(hexes, 4326))\nitems <-\n  stac(\"https://planetarycomputer.microsoft.com/api/stac/v1\") |>\n  stac_search(collections = \"io-biodiversity\",\n              bbox = c(box),\n              limit = 1000) |>\n  post_request() |>\n  items_sign(sign_fn = sign_planetary_computer())\n```\n:::\n\n\nCreated a desired cube view (Note, for simplicitly I show this at a much coarser resolution than the native data, allowing the gdalwarper to aggregate so this will run much faster)\n\n\n::: {.cell}\n\n```{.r .cell-code}\ncube <- gdalcubes::stac_image_collection(items$features, asset_names = \"data\")\nv <- cube_view(srs = \"EPSG:4326\",\n               extent = list(t0 = \"2020-01-01\", t1 = \"2020-12-31\",\n                             left = box[1], right = box[3],\n                             top = box[4], bottom = box[2]),\n               nx = 1000, ny=1000, dt = \"P1Y\")\n\nQ <- raster_cube(cube,v)\n```\n:::\n\n\nQuick plot of the data at requested cube resolution.  \n\n\n::: {.cell}\n\n```{.r .cell-code}\nQ |> plot(col=viridisLite::viridis(10))\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-7-1.png){width=672}\n:::\n:::\n\n\nExtract mean value of the BII for each hex.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nbii_hexes <- Q |> gdalcubes::extract_geom(hexes, FUN=mean)\n```\n:::\n\n\n\nPlot\n\n\n::: {.cell}\n\n```{.r .cell-code}\nbii_hexes <- hexes |>\n  tibble::rowid_to_column(\"FID\") |> \n  left_join(bii_hexes)\n\ntm_shape(bii_hexes) + tm_polygons(fill = \"data\", col=NA)\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-9-1.png){width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\ntmap_mode(\"plot\")\nsf = st_bbox(c(xmin=-123, ymin=39, xmax=-122, ymax=37.5), crs=4326)\ntm_shape(bii_hexes, bbox = sf) + tm_polygons(fill = \"data\", col=NA)\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-10-1.png){width=672}\n:::\n:::\n\n\n\n## Leaflet map\n\n\n::: {.cell}\n\n```{.r .cell-code}\ntmap_mode(\"view\")\ntm_shape(bii_hexes) + tm_polygons(fill = \"data\", col=NA)\n```\n:::\n",
+    "markdown": "---\ntitle: ACE Data\n---\n\n\n\nCalifornia Department of Fish and Wildlife (CDFW) provides a series of data products to identify \"areas of conservation emphasis\" (ACE), which are currently being used as part of [California's 30x30 conservation initiative](https://www.californianature.ca.gov/pages/30x30). At a recent event at the California Academy, one of the CDFW scientists mentioned it would be nice to see how the Biodiversity Intactness Index (BII) values compare to the CDFW prioritization layers. The BII (sometimes called the local or LBII) is one of a handful of indicators of global biodiversity change highlighted by [GEO BON](https://geobon.org/ebvs/indicators/). Various attempts to estimate this indicator have been made over the past decades -- the most recent and most high-resolution one currently available (and updated annually) is the Impact Observatory product based on the approach of Newbold et al and the PREDICTS database, which is [available in cloud-optimized format](https://planetarycomputer.microsoft.com/dataset/io-biodiversity) from the Planetary Computer STAC catalog. \n\nAside from drowning in acronyms, this task provides a great chance to apply the same tools seen in our [intro](intro.html) at a larger spatial scale.  We will average the rasterized 100m projections from the BII over each of the 63,890 polygon tiles used in the CDFW data. Even though we are now dealing with a raster layer that involves a very derived product (BII) from a different provider (Planetary Computer), and a much larger set of polygons from a different source and in different scale and projection, the process is almost identical to the intro example.  \n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(stars)\nlibrary(rstac)\nlibrary(gdalcubes)\nlibrary(tmap)\nlibrary(dplyr)\n\ngdalcubes::gdalcubes_options(parallel=24*2)\n```\n:::\n\n\nCDFW uses a hex-based tiling scheme for their data products. \nNote that these pre-date the modern [H3 hex tiling system](https://h3geo.org/), so do not provide the magic of [hierarchical (zoom level) tiles](https://h3geo.org/docs/highlights/indexing).\nCDFW makes the [ACE GIS Data](https://wildlife.ca.gov/Data/Analysis/Ace) freely available, but in one large zip archive. \nThe ACE data includes may different layers all using the same hex-based tiles. \nThe layer we pull here draws from their summary rankings on Terrestrial Irreplacability of species biodiversity.\n\nHere I'll use a public mirror on a Berkeley-based server for a more cloud-native access pattern so we don't have to download the whole thing.  We also index this in a simple [stac catalog](https://radiantearth.github.io/stac-browser/#/external/raw.githubusercontent.com/schmidtDSE/biodiversity-catalog/main/stac/v1/collections/cdfw_ace/summary.json?.asset=asset-ds2715)\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nurl <- \"/vsicurl/https://minio.carlboettiger.info/public-biodiversity/ACE_summary/ds2715.gdb\"\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nace <- st_read(url) \n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\ntm_shape(ace) + tm_polygons(fill = \"RwiRankEco\", col=NA)\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-5-1.png){width=672}\n:::\n:::\n\n\nFind Biodiversity Intactness Index COG tiles from Planetary Computer using STAC search:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nbox <- st_bbox(st_transform(ace, 4326))\nitems <-\n  stac(\"https://planetarycomputer.microsoft.com/api/stac/v1\") |>\n  stac_search(collections = \"io-biodiversity\",\n              bbox = c(box),\n              limit = 1000) |>\n  post_request() |>\n  items_sign(sign_fn = sign_planetary_computer())\n```\n:::\n\n\nCreated a desired cube view \n\n\n::: {.cell}\n\n```{.r .cell-code}\ncube <- gdalcubes::stac_image_collection(items$features, asset_names = \"data\")\nbox <- st_bbox(ace)\nv <- cube_view(srs = \"EPSG:3310\",\n               extent = list(t0 = \"2020-07-01\", t1 = \"2020-07-31\",\n                                left = box[1], right = box[3],\n                                top = box[4], bottom = box[2]),\n               dx = 100, dy=100, dt = \"P1M\")\n\nQ <- raster_cube(cube,v)\n```\n:::\n\n\nQuick plot of the data at requested cube resolution.  \n\n\n::: {.cell}\n\n```{.r .cell-code}\nQ |> plot(col=viridisLite::viridis(30), nbreaks=31)\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-8-1.png){width=672}\n:::\n:::\n\n\nExtract mean value of the BII for each hex.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nbii <- Q |> gdalcubes::extract_geom(ace, FUN=mean)\n```\n:::\n\n\n\n## Plots\n\nView the average BII index values by ACE hexagon:\n\n\n::: {.cell}\n\n```{.r .cell-code}\ntmap_mode(\"plot\")\n\n# color palette\nviridis <- tm_scale_continuous(values = viridisLite::viridis(30))\n\nbii <- ace |>\n  tibble::rowid_to_column(\"FID\") |> \n  left_join(bii)\n\ntm_shape(bii) + \n  tm_polygons(fill = \"data\", col=NA, fill.scale = viridis)\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-10-1.png){width=672}\n:::\n:::\n\n\nLet's zoom in manually:\n\n\n::: {.cell}\n\n```{.r .cell-code}\ntmap_mode(\"plot\")\nsf = st_bbox(c(xmin=-123, ymin=38.5, xmax=-122, ymax=37.5), crs=4326)\ntm_shape(bii, bbox = sf) + \n  tm_polygons(fill = \"data\", col=NA,\n              fill.scale = viridis)\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-11-1.png){width=672}\n:::\n:::\n\n\nNow that we have the data in convenient and familiar data structures, it is easy to analyze.\nOn average, a hexagon's irreplacability rank shows little correlation with BII:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nbii |> as_tibble() |> select(-Shape) |>\n  group_by(RwiRankEco) |> \n  summarise(BII = mean(data, na.rm=TRUE),\n            sd = sd(data, na.rm=TRUE)) |>\n  ggplot(aes(RwiRankEco, BII)) + \n  geom_col(aes(fill=RwiRankEco)) + \n  geom_linerange(aes(ymin = BII-2*sd, ymax = BII+2*sd), col = \"grey50\")\n```\n\n::: {.cell-output-display}\n![](ace_files/figure-html/unnamed-chunk-12-1.png){width=672}\n:::\n:::\n\n\n\n## Leaflet map\n\nWe can also render an interactive leaflet plot where we can zoom in and out and toggle basemaps:\n\n\n\n```{.r .cell-code}\ntmap_mode(\"view\")\n\nmap <- \n  tm_shape(bii) +\n  tm_polygons(fill = \"data\", col=NA, fill.scale = viridis) +\n  tm_shape(ace) +\n  tm_polygons(fill = \"RwiRankEco\", col=NA)\ntmap::tmap_save(map, \"bii_hexes.html\")\n```\n\n\n\n[interactive map](bii_hexes.html)\n\n",
     "supporting": [
       "ace_files"
     ],

diff --git a/_freeze/contents/ace/figure-html/unnamed-chunk-10-1.png b/_freeze/contents/ace/figure-html/unnamed-chunk-10-1.png
diff --git a/_freeze/contents/ace/figure-html/unnamed-chunk-11-1.png b/_freeze/contents/ace/figure-html/unnamed-chunk-11-1.png
diff --git a/_freeze/contents/ace/figure-html/unnamed-chunk-12-1.png b/_freeze/contents/ace/figure-html/unnamed-chunk-12-1.png
diff --git a/...nts/ace/figure-html/unnamed-chunk-4-1.png → ...nts/ace/figure-html/unnamed-chunk-5-1.png b/...nts/ace/figure-html/unnamed-chunk-4-1.png → ...nts/ace/figure-html/unnamed-chunk-5-1.png
diff --git a/_freeze/contents/ace/figure-html/unnamed-chunk-7-1.png b/_freeze/contents/ace/figure-html/unnamed-chunk-7-1.png
diff --git a/_freeze/contents/ace/figure-html/unnamed-chunk-8-1.png b/_freeze/contents/ace/figure-html/unnamed-chunk-8-1.png
diff --git a/_freeze/contents/ace/figure-html/unnamed-chunk-9-1.png b/_freeze/contents/ace/figure-html/unnamed-chunk-9-1.png
diff --git a/_quarto.yml b/_quarto.yml
@@ -38,6 +38,9 @@ website:
           - text: Introduction
             icon: play-btn
             href: contents/intro.qmd
+          - text: Biodiversity Priorities
+            icon: tree-fill
+            href: contents/ace.qmd  
           - text: NASA EarthData
             icon: rocket
             href: contents/earthdata.qmd  

diff --git a/contents/ace.qmd b/contents/ace.qmd
@@ -1,7 +1,12 @@
 ---
-
+title: ACE Data
 ---
 
+
+California Department of Fish and Wildlife (CDFW) provides a series of data products to identify "areas of conservation emphasis" (ACE), which are currently being used as part of [California's 30x30 conservation initiative](https://www.californianature.ca.gov/pages/30x30). At a recent event at the California Academy, one of the CDFW scientists mentioned it would be nice to see how the Biodiversity Intactness Index (BII) values compare to the CDFW prioritization layers. The BII (sometimes called the local or LBII) is one of a handful of indicators of global biodiversity change highlighted by [GEO BON](https://geobon.org/ebvs/indicators/). Various attempts to estimate this indicator have been made over the past decades -- the most recent and most high-resolution one currently available (and updated annually) is the Impact Observatory product based on the approach of Newbold et al and the PREDICTS database, which is [available in cloud-optimized format](https://planetarycomputer.microsoft.com/dataset/io-biodiversity) from the Planetary Computer STAC catalog. 
+
+Aside from drowning in acronyms, this task provides a great chance to apply the same tools seen in our [intro](intro.html) at a larger spatial scale.  We will average the rasterized 100m projections from the BII over each of the 63,890 polygon tiles used in the CDFW data. Even though we are now dealing with a raster layer that involves a very derived product (BII) from a different provider (Planetary Computer), and a much larger set of polygons from a different source and in different scale and projection, the process is almost identical to the intro example.  
+
 ```{r include=FALSE}
 knitr::opts_chunk$set(message = FALSE, warning = FALSE)
 ```
@@ -12,12 +17,15 @@ library(rstac)
 library(gdalcubes)
 library(tmap)
 library(dplyr)
-tmap_mode("plot")
 
 gdalcubes::gdalcubes_options(parallel=24*2)
 ```
 
-CDFW makes the [ACE GIS Data](https://wildlife.ca.gov/Data/Analysis/Ace) freely available, but in one large zip archive.
+CDFW uses a hex-based tiling scheme for their data products. 
+Note that these pre-date the modern [H3 hex tiling system](https://h3geo.org/), so do not provide the magic of [hierarchical (zoom level) tiles](https://h3geo.org/docs/highlights/indexing).
+CDFW makes the [ACE GIS Data](https://wildlife.ca.gov/Data/Analysis/Ace) freely available, but in one large zip archive. 
+The ACE data includes may different layers all using the same hex-based tiles. 
+The layer we pull here draws from their summary rankings on Terrestrial Irreplacability of species biodiversity.
 
 Here I'll use a public mirror on a Berkeley-based server for a more cloud-native access pattern so we don't have to download the whole thing.  We also index this in a simple [stac catalog](https://radiantearth.github.io/stac-browser/#/external/raw.githubusercontent.com/schmidtDSE/biodiversity-catalog/main/stac/v1/collections/cdfw_ace/summary.json?.asset=asset-ds2715)
 
@@ -26,17 +34,20 @@ Here I'll use a public mirror on a Berkeley-based server for a more cloud-native
 url <- "/vsicurl/https://minio.carlboettiger.info/public-biodiversity/ACE_summary/ds2715.gdb"
 ```
 
-CDFW uses hexes!  Note that these pre-date the modern [H3 hex tiling system](https://h3geo.org/), so do not provide the magic of [hierarchical (zoom level) tiles](https://h3geo.org/docs/highlights/indexing)
+
+
+```{r results="hide"}
+ace <- st_read(url) 
+```
 
 ```{r}
-hexes <- st_read(url) 
-tm_shape(hexes) + tm_polygons(fill = "RwiRankEco", col=NA)
+tm_shape(ace) + tm_polygons(fill = "RwiRankEco", col=NA)
 ```
 
 Find Biodiversity Intactness Index COG tiles from Planetary Computer using STAC search:
 
 ```{r}
-box <- st_bbox(st_transform(hexes, 4326))
+box <- st_bbox(st_transform(ace, 4326))
 items <-
   stac("https://planetarycomputer.microsoft.com/api/stac/v1") |>
   stac_search(collections = "io-biodiversity",
@@ -47,54 +58,93 @@ items <-
 
 ```
 
-Created a desired cube view (Note, for simplicitly I show this at a much coarser resolution than the native data, allowing the gdalwarper to aggregate so this will run much faster)
+Created a desired cube view 
 
 ```{r}
 cube <- gdalcubes::stac_image_collection(items$features, asset_names = "data")
-v <- cube_view(srs = "EPSG:4326",
-               extent = list(t0 = "2020-01-01", t1 = "2020-12-31",
-                             left = box[1], right = box[3],
-                             top = box[4], bottom = box[2]),
-               nx = 1000, ny=1000, dt = "P1Y")
+box <- st_bbox(ace)
+v <- cube_view(srs = "EPSG:3310",
+               extent = list(t0 = "2020-07-01", t1 = "2020-07-31",
+                                left = box[1], right = box[3],
+                                top = box[4], bottom = box[2]),
+               dx = 100, dy=100, dt = "P1M")
 
 Q <- raster_cube(cube,v)
 ```
 
 Quick plot of the data at requested cube resolution.  
 
 ```{r}
-Q |> plot(col=viridisLite::viridis(10))
+Q |> plot(col=viridisLite::viridis(30), nbreaks=31)
 ```
 
 Extract mean value of the BII for each hex.
 
 ```{r}
-bii_hexes <- Q |> gdalcubes::extract_geom(hexes, FUN=mean)
+bii <- Q |> gdalcubes::extract_geom(ace, FUN=mean)
 ```
 
 
-Plot
+## Plots
+
+View the average BII index values by ACE hexagon:
 
 ```{r}
-bii_hexes <- hexes |>
+tmap_mode("plot")
+
+# color palette
+viridis <- tm_scale_continuous(values = viridisLite::viridis(30))
+
+bii <- ace |>
   tibble::rowid_to_column("FID") |> 
-  left_join(bii_hexes)
+  left_join(bii)
 
-tm_shape(bii_hexes) + tm_polygons(fill = "data", col=NA)
+tm_shape(bii) + 
+  tm_polygons(fill = "data", col=NA, fill.scale = viridis)
 ```
 
+Let's zoom in manually:
+
 ```{r}
 tmap_mode("plot")
-sf = st_bbox(c(xmin=-123, ymin=39, xmax=-122, ymax=37.5), crs=4326)
-tm_shape(bii_hexes, bbox = sf) + tm_polygons(fill = "data", col=NA)
+sf = st_bbox(c(xmin=-123, ymin=38.5, xmax=-122, ymax=37.5), crs=4326)
+tm_shape(bii, bbox = sf) + 
+  tm_polygons(fill = "data", col=NA,
+              fill.scale = viridis)
+
+```
+
+Now that we have the data in convenient and familiar data structures, it is easy to analyze.
+On average, a hexagon's irreplacability rank shows little correlation with BII:
 
+```{r}
+library(tidyverse)
+bii |> as_tibble() |> select(-Shape) |>
+  group_by(RwiRankEco) |> 
+  summarise(BII = mean(data, na.rm=TRUE),
+            sd = sd(data, na.rm=TRUE)) |>
+  ggplot(aes(RwiRankEco, BII)) + 
+  geom_col(aes(fill=RwiRankEco)) + 
+  geom_linerange(aes(ymin = BII-2*sd, ymax = BII+2*sd), col = "grey50")
+  
 ```
 
 
 ## Leaflet map
 
-```{r}
+We can also render an interactive leaflet plot where we can zoom in and out and toggle basemaps:
+
+```{r results = "asis"}
 tmap_mode("view")
-tm_shape(bii_hexes) + tm_polygons(fill = "data", col=NA)
+
+map <- 
+  tm_shape(bii) +
+  tm_polygons(fill = "data", col=NA, fill.scale = viridis) +
+  tm_shape(ace) +
+  tm_polygons(fill = "RwiRankEco", col=NA)
+tmap::tmap_save(map, "bii_hexes.html")
 ```
 
+
+[interactive map](bii_hexes.html)
+