Created pipeline to get zonal statistics from layers in the STAC catalog

pipelines/Zonal_statistics.json

@@ -0,0 +1,168 @@
+{
+  "nodes": [
+    {
+      "id": "0",
+      "type": "io",
+      "position": {
+        "x": 184.859375,
+        "y": 226
+      },
+      "data": {
+        "descriptionFile": "data>getBoundingBox.yml"
+      }
+    },
+    {
+      "id": "1",
+      "type": "io",
+      "position": {
+        "x": 815.859375,
+        "y": 262
+      },
+      "data": {
+        "descriptionFile": "data>zonalStats.yml"
+      }
+    },
+    {
+      "id": "2",
+      "type": "output",
+      "position": {
+        "x": 648.859375,
+        "y": 224
+      },
+      "data": {
+        "label": "Output"
+      }
+    },
+    {
+      "id": "3",
+      "type": "output",
+      "position": {
+        "x": 1221.859375,
+        "y": 263
+      },
+      "data": {
+        "label": "Output"
+      }
+    },
+    {
+      "id": "4",
+      "type": "output",
+      "position": {
+        "x": 1205.859375,
+        "y": 276
+      },
+      "data": {
+        "label": "Output"
+      }
+    }
+  ],
+  "edges": [
+    {
+      "source": "0",
+      "sourceHandle": "study_area_polygon",
+      "target": "2",
+      "targetHandle": null,
+      "id": "reactflow__edge-0study_area_polygon-2"
+    },
+    {
+      "source": "0",
+      "sourceHandle": "bbox",
+      "target": "1",
+      "targetHandle": "bbox",
+      "id": "reactflow__edge-0bbox-1bbox"
+    },
+    {
+      "source": "1",
+      "sourceHandle": "rasters",
+      "target": "3",
+      "targetHandle": null,
+      "id": "reactflow__edge-1rasters-3"
+    },
+    {
+      "source": "1",
+      "sourceHandle": "stats",
+      "target": "4",
+      "targetHandle": null,
+      "id": "reactflow__edge-1stats-4"
+    },
+    {
+      "source": "0",
+      "sourceHandle": "study_area_polygon",
+      "target": "1",
+      "targetHandle": "study_area_polygon",
+      "id": "reactflow__edge-0study_area_polygon-1study_area_polygon"
+    }
+  ],
+  "inputs": {
+    "data>getBoundingBox.yml@0|country": {
+      "description": "Country of interest",
+      "label": "Country",
+      "type": "text",
+      "weight": 0
+    },
+    "data>getBoundingBox.yml@0|state": {
+      "description": "State or province of interest",
+      "label": "State/Province",
+      "type": "text",
+      "weight": 1
+    },
+    "data>getBoundingBox.yml@0|study_area_file": {
+      "description": "Upload a file of a custom study area",
+      "label": "File for study area",
+      "type": "text",
+      "weight": 2
+    },
+    "data>zonalStats.yml@1|stac_url": {
+      "description": "URL of the Stac catalog",
+      "label": "stac_url",
+      "type": "text",
+      "example": "https://stac.geobon.org/",
+      "weight": 3
+    },
+    "data>zonalStats.yml@1|collections_items": {
+      "description": "Vector of strings, collection name followed by '|' followed by item id",
+      "label": "collections_items",
+      "type": "text[]",
+      "example": [
+        "bii_nhm|bii_nhm_10km_2000",
+        "bii_nhm|bii_nhm_10km_2005",
+        "bii_nhm|bii_nhm_10km_2010",
+        "bii_nhm|bii_nhm_10km_2015",
+        "bii_nhm|bii_nhm_10km_2020"
+      ],
+      "weight": 4
+    },
+    "data>zonalStats.yml@1|summary_statistic": {
+      "description": "Summary statistic for layers",
+      "label": "Summary statistic",
+      "type": "options",
+      "options": [
+        "mean",
+        "median",
+        "mode"
+      ],
+      "example": "mean",
+      "weight": 5
+    }
+  },
+  "outputs": {
+    "data>getBoundingBox.yml@0|study_area_polygon": {
+      "description": "Represents the map of the study area",
+      "label": "Polygon of study area",
+      "type": "application/geo+json",
+      "weight": 0
+    },
+    "data>zonalStats.yml@1|rasters": {
+      "description": "array of output raster paths",
+      "label": "Rasters",
+      "type": "image/tiff;application=geotiff[]",
+      "weight": 1
+    },
+    "data>zonalStats.yml@1|stats": {
+      "description": "Summary statistic in the polygon",
+      "label": "Summary statistic",
+      "type": "text/csv",
+      "weight": 2
+    }
+  }
+}

scripts/data/getBoundingBox.R

@@ -0,0 +1,41 @@
+# Script to get the bounding box from a country or region
+
+library("rjson")
+library("sf")
+
+input <- fromJSON(file=file.path(outputFolder, "input.json"))
+
+# define study area
+if (is.null(input$studyarea_file)){ # if there is no study area file input
+  if (is.null(input$state)){ # if there is only a country input (no state) 
+    input$country <- gsub(" ", "+", input$country) # Change spaces to + signs to work in the URL 
+    study_area<- paste0("https://geoio.biodiversite-quebec.ca/country_geojson/?country_name=", input$country) # study area url 
+  } else { # if a state is defined
+   input$country <- gsub(" ", "+", input$country)
+   input$state <- gsub(" ", "+", input$state)
+    study_area<- paste0("https://geoio.biodiversite-quebec.ca/state_geojson/?country_name=", input$country, "&state_name=", input$state)
+  } } else {study_area <- input$studyarea_file}
+
+study_area_polygon<- sf::st_read(study_area)  # load study area as sf object
+
+if(nrow(study_area_polygon)==0){
+  stop("Study area polygon does not exist. Check spelling of country and state names. Check if region contains protected areas")
+}  # stop if object is empty
+
+# Save study area and protected area data
+study_area_polygon_path<- file.path(outputFolder, "study_area_polygon.geojson") # Define the file path for the protected area polygon output
+sf::st_write(study_area_polygon, study_area_polygon_path, delete_dsn = T)
+
+
+# create bounding box
+bbox <- sf::st_bbox(study_area_polygon)
+
+
+bbox<-unname(st_bbox(study_area_polygon))
+
+output <- list("bbox" = bbox, "study_area_polygon" = study_area_polygon_path)
+print(output)
+
+### return outpu
+jsonData <- toJSON(output, indent=2)
+write(jsonData, file.path(outputFolder,"output.json"))

scripts/data/getBoundingBox.yml

@@ -0,0 +1,28 @@
+script: getBoundingBox.R
+name: Get bounding box
+description: "Extract items and calculate zonal statistics from various collections on the GEO BON Stac Catalog."
+author:
+  - name: Jory Griffith
+    identifier: https://orcid.org/0000-0001-6020-6690 
+inputs:
+  country:
+    label: Country
+    description: Country of interest
+    type: text
+  state:
+    label: State/Province
+    description: State or province of interest
+    type: text
+  study_area_file:
+    label: File for study area
+    description: Upload a file of a custom study area
+    type: text
+outputs:
+  study_area_polygon:
+    label: Polygon of study area
+    description: Represents the map of the study area
+    type: application/geo+json
+  bbox:
+    label: bbox
+    description: boundary box around area of interest
+    type: float[]

scripts/data/loadFromStac.yml

@@ -60,6 +60,7 @@ conda:
   dependencies:
     - libgdal
     - proj
+    - r-rjson
     - r-proj
     - r-gdalcubes
     - r-rstac

scripts/data/loadFromStacFun.R

@@ -82,6 +82,7 @@ load_cube <-
       feats<-it_obj$features
     }
     print(ids)
+    print(feats)
     if (!is.null(variable)) {
       print("Variable is null")
       st <- gdalcubes::stac_image_collection(

scripts/data/zonalStats.R

@@ -12,42 +12,67 @@ source(paste(Sys.getenv("SCRIPT_LOCATION"), "/data/loadFromStacFun.R", sep = "/"
 
 input <- fromJSON(file=file.path(outputFolder, "input.json"))
 
-# create bounding box
-bbox <- sf::st_bbox(c(xmin = input$bbox[1], ymin = input$bbox[2],
-            xmax = input$bbox[3], ymax = input$bbox[4]), crs = sf::st_crs(input$proj)) 
-
+output<- tryCatch({
 # Collections items
 if (length(input$collections_items) == 0) {
     stop('Please specify collections_items') # if no collections items are specified
 } else {
     collections_items <- input$collections_items
 }
 
-# define study area
-if (is.null(input$studyarea_file)){ # if there is no study area file input
-  if (is.null(input$state)){ # if there is only a country input (no state) 
-    input$country <- gsub(" ", "+", input$country) # Change spaces to + signs to work in the URL 
-    study_area<- paste0("https://geoio.biodiversite-quebec.ca/country_geojson/?country_name=", input$country) # study area url 
-  } else { # if a state is defined
-   input$country <- gsub(" ", "+", input$country)
-   input$state <- gsub(" ", "+", input$state)
-    study_area<- paste0("https://geoio.biodiversite-quebec.ca/state_geojson/?country_name=", input$country, "&state_name=", input$state)
-  } } else {study_area <- input$studyarea_file}
-
-# Read in study area polygon
-study_area_polygon<- sf::st_read(study_area)  # load study area as sf object
+# Load study area
+study_area <- st_read(input$study_area_polygon) %>% st_transform("EPSG:4326")
 
+# Define bounding box
+bbox <- sf::st_bbox(c(xmin = input$bbox[1], ymin = input$bbox[2],
+            xmax = input$bbox[3], ymax = input$bbox[4]), crs=st_crs(4326)) 
 
 # Load cube using the loadFromStacFun
-dat_cube <- load_cube(stac_path=input$stac_url, collections=collections_items, bbox=bbox,
-                        srs.cube = input$proj)
+raster_layers <- list()
+stats_list <- list()
+for(i in 1:length(collections_items)){
 
+  ci <- strsplit(collections_items[i], split = "|", fixed=TRUE)[[1]] # split into collection and layers
+
+  dat_cube <- load_cube(stac_path=input$stac_url, collections=ci[1], ids=ci[2], bbox=bbox,
+                           srs.cube = "EPSG:4326", layers=NULL, variable=NULL,
+                            t0 = NULL,   t1 = NULL, temporal.res = "P1D", aggregation = "mean",  resampling = "near")
+raster_layers[[i]]=dat_cube # add to a list of raster layers
+names(raster_layers)[i] <- ci[2]
+
+if(input$summary_statistic == "mean"){
+stats_each <- gdalcubes::extract_geom(cube=dat_cube, sf=study_area, FUN=mean, reduce_time=TRUE)
+} else if(input$summary_statistic == "median"){
+  stats_each <- gdalcubes::extract_geom(cube=dat_cube, sf=study_area, FUN=median, reduce_time=TRUE)
+} else {
+  stats_each <- gdalcubes::extract_geom(cube=dat_cube, sf=study_area, FUN=mode, reduce_time=TRUE)}
+
+stats_each$name <- ci[2]
+stats_list[[i]] = as.data.frame(stats_each)
+}
+
+stats <- bind_rows(stats_list)
+stats <- stats[,c(3,2)]
+names(stats)[names(stats) == 'data']<-input$summary_statistic
 
-# Calculate summary statistics - mean
-stats <- gdalcubes::extract_geom(cube=dat_cube, sf=study_area_polygon, FUN=mean)
 stats_path <- file.path(outputFolder, "stats.csv")
+write.csv(stats, stats_path, row.names=F)
+
+layer_paths<-c()
+# need to rename these to be more descriptive
+output_raster_layers <- file.path(outputFolder)
+
+for (i in 1:length(raster_layers)) {
+  ff <- tempfile(pattern = paste0(names(raster_layers[i]),'_'))
+  out<-gdalcubes::write_tif(raster_layers[i][[1]], dir = output_raster_layers, prefix=basename(ff), creation_options = list("COMPRESS" = "DEFLATE"), COG=TRUE, write_json_descr=TRUE)
+  fp <- paste0(out[1])
+  layer_paths <- cbind(layer_paths,fp)
+}
+
+
+output <- list("stats" = stats_path, "rasters" = layer_paths)
+}, error = function(e) { list(error = conditionMessage(e)) })
 
-output <- list("stats" = stats_path)
 jsonData <- toJSON(output, indent=2)
 write(jsonData, file.path(outputFolder,"output.json"))