diff --git a/vignettes/04-synthesis-data.Rmd b/vignettes/04-synthesis-data.Rmd index d23529b..485f7ca 100644 --- a/vignettes/04-synthesis-data.Rmd +++ b/vignettes/04-synthesis-data.Rmd @@ -5,10 +5,6 @@ The objective of this vignette is to walk through how to combine our several typ For the first analysis, we want to figure out how the number of sufficiently warm days affects the amount of canopy cover at our site. We do this by combining the canopy cover data with the meteorological data on growing degree days, then modeling and plotting their relationship. We are specifically interested in figuring out when the increase in canopy cover starts to slow down in response to warm temperature days. The second analysis compares greenness from image data with canopy cover. -The second analysis uses the gdal_translate tool from the [GDAL](https://www.gdal.org/) package. -You can use one of the prepared downloads availble on the GDAL web site to install the needed software tool. -Alternatively, it's possible to download and build the tools on your system. -For example, on the Mac, the command is `brew install gdal`. ## Get and join data @@ -24,6 +20,8 @@ library(jsonlite) library(lubridate) library(traits) library(inflection) +library(sf) +library(stringr) options(betydb_url = "https://terraref.ncsa.illinois.edu/bety/", betydb_api_version = 'v1') ``` @@ -131,59 +129,28 @@ ggplot(data.frame(inf_points = unique(all_cultivars$inf_point))) + ## Get image data In this examnple we will extract our plot data from a series of images taken in May of Season 6, measure its "greeness" annd plot that against the plant heights from above in this vignette. + The chosen statistic here is the normalised green-red difference index, NGRDI=(R-G)/(R+G) (Rasmussen et al., 2016), which uses the red and green bands from the image raster. Below we retrieve all the available plots for a particular date, then find and convert the plot boundary JSON into tuples. We will use these tuples to extract the data for our plot. ```{r get_plot_boundary} -library(traits) -library(stringr) - -# Function for breaking apart a corner into its Lat, Lon components -getLatLon <- function(corner){ - p <- strsplit(corner, ' ') - return (c(p[[1]][1], p[[1]][2])) -} - -# Gets the bounding box of the array of points -getBounds <-function(bounds){ - minX <- NA - minY <- NA - maxX <- NA - maxY <- NA - for (c in unique(bounds)){ - p = getLatLon(c) - if (is.na(minX) || (minX > p[2])) - minX = p[2] - if (is.na(minY) || (minY > p[1])) - minY = p[1] - if (is.na(maxX) || (maxX < p[2])) - maxX = p[2] - if (is.na(maxY) || (maxY < p[1])) - maxY = p[1] - } - - return (c(minX, minY, maxX, maxY)) -} - # Setting up our options options(betydb_url = "https://terraref.ncsa.illinois.edu/bety/", betydb_api_version = 'v1') -# Makiong the query for our site +# Making the query for our site sites <- betydb_query(table = "sites", sitename = "MAC Field Scanner Season 6 Range 19 Column 1") # Assigning the geometry of the site (GeoJSON format) site.geom <- sites$geometry -# Stripping out the extra information to get to the points -complete_str <- str_match_all(site.geom, '(\\(\\(\\((.*)\\)\\)\\))')[[1]][, 3] -bounds <- strsplit(complete_str, ', ')[[1]] - -# Getting the bounding box of the polygon -bounding_box = getBounds(bounds) +# Convert the polygon to something we can clip with. CRS value represents WGS84 Lat/Long +site.shape <- st_as_sfc(site.geom,crs = 4326) +site.poly <- st_cast(site.shape, "POINT") +site.clip <- as(site.poly,"Spatial") ``` These are the names of the full field RGB data for the month of May. @@ -196,7 +163,6 @@ image_files <- c('fullfield_L1_ua-mac_2018-05-01_rgb_stereovis_ir_sensors_fullfield_sorghum6_shade_may2018_thumb.tif', 'fullfield_L1_ua-mac_2018-05-02_rgb_stereovis_ir_sensors_fullfield_sorghum6_shade_may2018_thumb.tif', 'fullfield_L1_ua-mac_2018-05-03_rgb_stereovis_ir_sensors_fullfield_sorghum6_shade_may2018_thumb.tif', - 'fullfield_L1_ua-mac_2018-05-04_rgb_stereovis_ir_sensors_fullfield_sorghum6_settingstest_may2018_thumb.tif', 'fullfield_L1_ua-mac_2018-05-05_rgb_stereovis_ir_sensors_fullfield_sorghum6_shade_may2018_thumb.tif', 'fullfield_L1_ua-mac_2018-05-06_rgb_stereovis_ir_sensors_fullfield_sorghum6_shade_may2018_thumb.tif', 'fullfield_L1_ua-mac_2018-05-08_rgb_stereovis_ir_sensors_fullfield_sorghum6_shade_may2018_thumb.tif', @@ -208,7 +174,6 @@ image_files <- 'fullfield_L1_ua-mac_2018-05-15_rgb_stereovis_ir_sensors_fullfield_sorghum6_sun_may2018_thumb.tif', 'fullfield_L1_ua-mac_2018-05-17_rgb_stereovis_ir_sensors_fullfield_sorghum6_shade_may2018_thumb.tif', 'fullfield_L1_ua-mac_2018-05-18_rgb_stereovis_ir_sensors_fullfield_sorghum6_sun_may2018_thumb.tif', - 'fullfield_L1_ua-mac_2018-05-19_rgb_stereovis_ir_sensors_plots_sorghum6_sun_thumb.tif', 'fullfield_L1_ua-mac_2018-05-20_rgb_stereovis_ir_sensors_plots_sorghum6_shade_thumb.tif', 'fullfield_L1_ua-mac_2018-05-21_rgb_stereovis_ir_sensors_fullfield_sorghum6_shade_may2018_thumb.tif', 'fullfield_L1_ua-mac_2018-05-22_rgb_stereovis_ir_sensors_plots_sorghum6_sun_thumb.tif', @@ -221,13 +186,11 @@ image_files <- image_files <- file.path("vignettes/", image_files) ``` - We will loop through these images, extract our plot data, and calculate the "greeness" of each extract. We are using the name of the file to extract the date for later. ```{r synth_get_greeness} library(raster) -library(stringr) # Extract the date from the file name getDate <- function(file_name){ @@ -235,39 +198,23 @@ getDate <- function(file_name){ return(date) } - -# Get the clip coordinates into the correct order -clip_coords <- paste(toString(bounding_box[4]),' ',toString(bounding_box[3]), - ' ',toString(bounding_box[2]),' ',toString(bounding_box[1])) - - -# Returns the greeness value of the clipped image +# Returns the greeness value of the plot in the specified file getGreeness <- function(file_name, clip_coords){ - out_file <- "extract.tif" - - # Execute the GDAL command to extract the plot - command = paste("gdal_translate -projwin ",clip_coords," ",file_name," ",out_file) - system(command) - - # Load the red & green bands of the image and calculate the greeness value - red_image <- raster(out_file, band = 1) - cellStats(red_image, stat = "mean") + + band_image <- raster(file_name, band = 1) + red_crop <- crop(band_image, clip_coords) - green_image <- raster(out_file, band = 2) - cellStats(green_image, stat = "mean") + band_image <- raster(file_name, band = 2) + green_crop <- crop(band_image, clip_coords) - add_rasters <- green_image + red_image + add_rasters <- green_crop + red_crop numerator <- cellStats(add_rasters, stat = "sum") - subtract_rasters <- green_image - red_image + subtract_rasters <- green_crop - red_crop denominator <- cellStats(subtract_rasters, stat = "sum") greeness <- numerator / denominator - # Remove the temporary file - if (file.exists(out_file)) - file.remove(out_file) - return(greeness) } @@ -275,7 +222,7 @@ getGreeness <- function(file_name, clip_coords){ day <- sapply(image_files, getDate, USE.NAMES = FALSE) # Extract all the greeness for the plot -greeness <- sapply(image_files, getGreeness, clip_coords=clip_coords, USE.NAMES = FALSE) +greeness <- sapply(image_files, getGreeness, clip_coords=site.clip, USE.NAMES = FALSE) # Build the final day and greeness greenness_df <- data.frame(day,greeness) @@ -292,7 +239,7 @@ trait_canopy_cover <- betydb_query(table = "search", trait = "canopy_cover", date = "~2018 May", limit = "none") - + trait_canopy_cover <- trait_canopy_cover %>% mutate(day = as.Date(raw_date)) ```