freeze 🧊

_freeze/contents/r-intro/execute-results/html.json

@@ -0,0 +1,16 @@
+{
+  "hash": "855ffee589b50d7f9a482b717bb12cb4",
+  "result": {
+    "markdown": "---\ntitle: \"Intro\"\nformat: html\n---\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\n```\n\n::: {.cell-output .cell-output-stderr}\n```\n── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──\n✔ dplyr     1.1.3     ✔ readr     2.1.4\n✔ forcats   1.0.0     ✔ stringr   1.5.0\n✔ ggplot2   3.4.3     ✔ tibble    3.2.1\n✔ lubridate 1.9.3     ✔ tidyr     1.3.0\n✔ purrr     1.0.2     \n── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──\n✖ dplyr::filter() masks stats::filter()\n✖ dplyr::lag()    masks stats::lag()\nℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors\n```\n:::\n\n```{.r .cell-code}\nlibrary(rstac)\nlibrary(gdalcubes)\nlibrary(stars)\n```\n\n::: {.cell-output .cell-output-stderr}\n```\nLoading required package: abind\nLoading required package: sf\nLinking to GEOS 3.12.0, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE\n```\n:::\n\n```{.r .cell-code}\nlibrary(tmap)\n```\n\n::: {.cell-output .cell-output-stderr}\n```\nThe legacy packages maptools, rgdal, and rgeos, underpinning the sp package,\nwhich was just loaded, will retire in October 2023.\nPlease refer to R-spatial evolution reports for details, especially\nhttps://r-spatial.org/r/2023/05/15/evolution4.html.\nIt may be desirable to make the sf package available;\npackage maintainers should consider adding sf to Suggests:.\nThe sp package is now running under evolution status 2\n     (status 2 uses the sf package in place of rgdal)\n\nAttaching package: 'tmap'\n\nThe following object is masked from 'package:datasets':\n\n    rivers\n```\n:::\n\n```{.r .cell-code}\ngdalcubes::gdalcubes_options(parallel = TRUE)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n## STAC Search over 400 million assets.\nbox <- c(xmin=-122.51006, ymin=37.70801, xmax=-122.36268, ymax=37.80668) \nstart_date <- \"2022-06-01\"\nend_date <- \"2022-08-01\"\nitems <- \n  stac(\"https://earth-search.aws.element84.com/v0/\") |>\n  stac_search(collections = \"sentinel-s2-l2a-cogs\",\n              bbox = box,\n              datetime = paste(start_date, end_date, sep=\"/\"),\n              limit = 100) |>\n  post_request() \n```\n:::\n\n\n## Low-level approach\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# 12 matching features\nlength(items$features)\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n[1] 12\n```\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nitems$features[[1]] |> names()\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n [1] \"type\"            \"stac_version\"    \"stac_extensions\" \"id\"             \n [5] \"bbox\"            \"geometry\"        \"properties\"      \"collection\"     \n [9] \"assets\"          \"links\"          \n```\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nitems$features[[1]]$assets |> names()\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n [1] \"thumbnail\" \"overview\"  \"info\"      \"metadata\"  \"visual\"    \"B01\"      \n [7] \"B02\"       \"B03\"       \"B04\"       \"B05\"       \"B06\"       \"B07\"      \n[13] \"B08\"       \"B8A\"       \"B09\"       \"B11\"       \"B12\"       \"AOT\"      \n[19] \"WVP\"       \"SCL\"      \n```\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nb04 <- items$features[[1]]$assets$B04 \nc(b04$title, b04$type)\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n[1] \"Band 4 (red)\"                                            \n[2] \"image/tiff; application=geotiff; profile=cloud-optimized\"\n```\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\ni <- 1\n\nb02 <- paste0(\"/vsicurl/\", items$features[[i]]$assets$B02$href)\nb03 <- paste0(\"/vsicurl/\", items$features[[i]]$assets$B03$href)\nb04 <- paste0(\"/vsicurl/\", items$features[[i]]$assets$B04$href)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nvis <- read_stars(c(b04,b03,b02))\nb04\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n[1] \"/vsicurl/https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2022/7/S2B_10SEG_20220728_0_L2A/B04.tif\"\n```\n:::\n\n```{.r .cell-code}\nplot(vis)\n```\n\n::: {.cell-output .cell-output-stderr}\n```\ndownsample set to 8\n```\n:::\n\n::: {.cell-output-display}\n![](r-intro_files/figure-html/unnamed-chunk-7-1.png){width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n# Don't try this, will definitely crash!\n#ggplot() + geom_stars(data = vis) + scale_fill_identity()\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nred <- read_stars(b04)\nb08 <- paste0(\"/vsicurl/\", items$features[[i]]$assets$B08$href)\nnir <- read_stars(b08)\nndvi <- (nir - red) / (nir + red)\nplot(ndvi)\n```\n\n::: {.cell-output .cell-output-stderr}\n```\ndownsample set to 8\n```\n:::\n\n::: {.cell-output-display}\n![](r-intro_files/figure-html/unnamed-chunk-9-1.png){width=672}\n:::\n:::\n\n\nThis approach Does not mask clouds and fill or average over other images -- we've used only one of the matching assets so far. While a single image here covers our area of interest (AOI),  in an expanded spatial analysis we would need to tile together overlapping images to ensure we cover the full AOI.  This also operates at native resolution, which can hurt performance when scaling to larger AOI.\n\n\n\n## Newer cloud-native approach\n\n\n::: {.cell}\n\n```{.r .cell-code}\ncol <-\n  stac_image_collection(items$features,\n                        asset_names = c(\"B02\", \"B03\", \"B04\",\"B08\", \"SCL\"),\n                        property_filter = \\(x) {x[[\"eo:cloud_cover\"]] < 20})\n```\n\n::: {.cell-output .cell-output-stderr}\n```\nWarning in stac_image_collection(items$features, asset_names = c(\"B02\", : STAC\nasset with name 'SCL' does not include eo:bands metadata and will be considered\nas a single band source\n```\n:::\n\n```{.r .cell-code}\ncube <- cube_view(srs = \"EPSG:4326\",  \n                  extent = list(t0 = start_date, t1 = end_date,\n                                left = box[1], right = box[3],\n                                top = box[4], bottom = box[2]),\n                  nx = 2400, ny = 2400, dt = \"P1M\",\n                  aggregation = \"median\", resampling = \"average\")\n\nS2.mask <- image_mask(\"SCL\", values=c(3,8,9)) # mask clouds and cloud shadows\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nraster_cube(col, cube, mask = S2.mask) |>\n  select_bands(c(\"B02\", \"B03\", \"B04\")) |>\n  aggregate_time(dt=\"P1M\") |>\n  plot(rgb=3:1)\n```\n\n::: {.cell-output-display}\n![](r-intro_files/figure-html/unnamed-chunk-11-1.png){width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nndvi <- raster_cube(col, cube, mask = S2.mask) |>\n  select_bands(c(\"B04\", \"B08\")) |>\n  apply_pixel(\"(B08-B04)/(B08+B04)\", \"NDVI\") |>\n  reduce_time(c(\"mean(NDVI)\")) |>\n  st_as_stars()\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\ntm_shape(ndvi) + \n  tm_raster(col.scale = tm_scale_continuous(values = viridisLite::mako(30)))\n```\n\n::: {.cell-output .cell-output-stderr}\n```\nstars object downsampled to 1200 by 1200 cells.\n```\n:::\n\n::: {.cell-output .cell-output-stderr}\n```\nWarning in value[[3L]](cond): could not rename the data.table\n```\n:::\n\n::: {.cell-output-display}\n![](r-intro_files/figure-html/unnamed-chunk-13-1.png){width=672}\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nCASanFrancisco1937 <- \"https://dsl.richmond.edu/panorama/redlining/static/downloads/geojson/CASanFrancisco1937.geojson\"\nredlines <- st_read(paste0(\"/vsicurl/\", CASanFrancisco1937)) |> st_make_valid()\n```\n\n::: {.cell-output .cell-output-stdout}\n```\nReading layer `CASanFrancisco1937' from data source \n  `/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/downloads/geojson/CASanFrancisco1937.geojson' \n  using driver `GeoJSON'\nSimple feature collection with 97 features and 4 fields\nGeometry type: MULTIPOLYGON\nDimension:     XY\nBounding box:  xmin: -122.5101 ymin: 37.70801 xmax: -122.3627 ymax: 37.80668\nGeodetic CRS:  WGS 84\n```\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\ntm_shape(ndvi) + \n  tm_raster(col.scale = tm_scale_continuous(values = viridisLite::mako(30))) +\n  tm_shape(redlines) + tm_borders(col=\"holc_grade\", lwd = 4, \n                                  col.scale = tm_scale_categorical(values=c(\n                                    A = \"#729765\",\n                                    B = \"#75a4b3\",\n                                    C = \"#c6c167\",\n                                    D = \"#b0707b\" ))) +\n  tm_text(text = \"holc_grade\", size = 0.5, col = \"holc_grade\")\n```\n\n::: {.cell-output .cell-output-stderr}\n```\nstars object downsampled to 1200 by 1200 cells.\n```\n:::\n\n::: {.cell-output .cell-output-stderr}\n```\nWarning in value[[3L]](cond): could not rename the data.table\n```\n:::\n\n::: {.cell-output .cell-output-stderr}\n```\nWarning: (Set of) legends is too high and are therefore rescaled.\n```\n:::\n\n::: {.cell-output-display}\n![](r-intro_files/figure-html/unnamed-chunk-15-1.png){width=672}\n:::\n:::\n\n::: {.cell}\n\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n# more scalable performance with gdalcubes::extract_geom\nndvi_aves <- raster_cube(col, cube, mask = S2.mask) |>\n  select_bands(c(\"B04\", \"B08\")) |>\n  apply_pixel(\"(B08-B04)/(B08+B04)\", \"NDVI\") |>\n  reduce_time(c(\"mean(NDVI)\")) |>\n  extract_geom(redlines, FUN=mean)\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\nave_ndvi <- redlines |> \n  rowid_to_column(\"FID\") |>\n  left_join(ndvi_aves) \n```\n\n::: {.cell-output .cell-output-stderr}\n```\nJoining with `by = join_by(FID)`\n```\n:::\n\n```{.r .cell-code}\nave_ndvi |>\n  as_tibble() |>\n  group_by(holc_grade) |>\n  summarise(mean = mean(NDVI_mean))\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n# A tibble: 4 × 2\n  holc_grade  mean\n  <chr>      <dbl>\n1 A          0.316\n2 B          0.212\n3 C          0.195\n4 D          0.194\n```\n:::\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\n# \"standard method\"\n# aves <- stars::st_extract(ndvi, redlines, FUN=mean)\n\n# vec <- as_tibble(aves) |> left_join(redlines)\n#vec |> group_by(holc_grade) |>\n#   summarise(ndvi = mean(NDVI_mean, na.rm=TRUE))\n```\n:::\n\n::: {.cell}\n\n```{.r .cell-code}\ntm_basemap(server =providers$OpenStreetMap, zoom=13) + \ntm_shape(ave_ndvi) + \n  tm_polygons(\"NDVI_mean\", \n              fill.scale = tm_scale_continuous(values = \"Greens\")) +\n  tm_shape(redlines) + \n  tm_borders(col=\"holc_grade\", \n             lwd = 4, \n             col.scale = tm_scale_categorical(values=c(\n                A = \"#729765\",\n                B = \"#75a4b3\",\n                C = \"#c6c167\",\n                D = \"#b0707b\" ))) +\n  tm_text(text = \"holc_grade\", size = 0.5)\n```\n\n::: {.cell-output .cell-output-stderr}\n```\nSpatRaster object downsampled to 889 by 1125 cells.\n```\n:::\n\n::: {.cell-output-display}\n![](r-intro_files/figure-html/unnamed-chunk-20-1.png){width=672}\n:::\n:::",
+    "supporting": [
+      "r-intro_files"
+    ],
+    "filters": [
+      "rmarkdown/pagebreak.lua"
+    ],
+    "includes": {},
+    "engineDependencies": {},
+    "preserve": {},
+    "postProcess": true
+  }
+}