From 5c432d0e8f81ec8ba98c10afe67125df5600836c Mon Sep 17 00:00:00 2001 From: Chris Schnaufer Date: Fri, 8 Feb 2019 11:21:42 -0700 Subject: [PATCH 1/3] Removed gdal_transform command line tool requirement --- vignettes/04-synthesis-data.Rmd | 101 ++++++++------------------------ 1 file changed, 24 insertions(+), 77 deletions(-) diff --git a/vignettes/04-synthesis-data.Rmd b/vignettes/04-synthesis-data.Rmd index d23529b..8afc3a4 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 @@ -18,12 +14,15 @@ The second is the cumulative growing degree days for all of 2017, which were cal They are combined by their common column, the date. ```{r synth_setup} +library(raster) library(dplyr) library(ggplot2) 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') ``` @@ -137,53 +136,22 @@ Below we retrieve all the available plots for a particular date, then find and c 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') + betydb_api_version = 'v1', + betydb_key = readLines('.betykey')) -# 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. @@ -217,57 +185,37 @@ image_files <- ) ``` -```{r, echo = F} -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){ date <- str_match_all(file_name, '[0-9]{4}-[0-9]{2}-[0-9]{2}')[[1]][,1] 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) + band_image = raster(file_name, band = 1) + red_crop = crop(band_image, clip_coords) - # 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 = 2) + green_crop = crop(band_image, clip_coords) - green_image <- raster(out_file, band = 2) - cellStats(green_image, stat = "mean") - - 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) - + + # If we can't get a greeness value we return zero + if (is.nan(greeness)) { + greeness = 0; + } + return(greeness) } @@ -275,7 +223,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 +240,6 @@ 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)) ``` From 10ef1ef6360306bb54f4f4686c2519e48c1cce9d Mon Sep 17 00:00:00 2001 From: Chris Schnaufer Date: Fri, 8 Feb 2019 11:39:54 -0700 Subject: [PATCH 2/3] Fixing an incomplete manual merge --- vignettes/04-synthesis-data.Rmd | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/vignettes/04-synthesis-data.Rmd b/vignettes/04-synthesis-data.Rmd index 8afc3a4..42d5fe6 100644 --- a/vignettes/04-synthesis-data.Rmd +++ b/vignettes/04-synthesis-data.Rmd @@ -138,8 +138,7 @@ We will use these tuples to extract the data for our plot. ```{r get_plot_boundary} # Setting up our options options(betydb_url = "https://terraref.ncsa.illinois.edu/bety/", - betydb_api_version = 'v1', - betydb_key = readLines('.betykey')) + betydb_api_version = 'v1') # Making the query for our site sites <- betydb_query(table = "sites", @@ -185,6 +184,10 @@ image_files <- ) ``` +```{r, echo = F} +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. @@ -240,6 +243,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)) ``` From e4f3134d86ebed3c0bd856ec7e6db1899e60730f Mon Sep 17 00:00:00 2001 From: Chris Schnaufer Date: Fri, 8 Feb 2019 12:55:49 -0700 Subject: [PATCH 3/3] Changed from pull request --- vignettes/04-synthesis-data.Rmd | 28 ++++++++++++---------------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/vignettes/04-synthesis-data.Rmd b/vignettes/04-synthesis-data.Rmd index 42d5fe6..485f7ca 100644 --- a/vignettes/04-synthesis-data.Rmd +++ b/vignettes/04-synthesis-data.Rmd @@ -14,7 +14,6 @@ The second is the cumulative growing degree days for all of 2017, which were cal They are combined by their common column, the date. ```{r synth_setup} -library(raster) library(dplyr) library(ggplot2) library(jsonlite) @@ -130,6 +129,7 @@ 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. @@ -148,9 +148,9 @@ sites <- betydb_query(table = "sites", site.geom <- sites$geometry # 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") +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. @@ -163,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', @@ -175,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', @@ -192,6 +190,8 @@ We will loop through these images, extract our plot data, and calculate the "gre We are using the name of the file to extract the date for later. ```{r synth_get_greeness} +library(raster) + # Extract the date from the file name getDate <- function(file_name){ date <- str_match_all(file_name, '[0-9]{4}-[0-9]{2}-[0-9]{2}')[[1]][,1] @@ -200,11 +200,12 @@ getDate <- function(file_name){ # Returns the greeness value of the plot in the specified file getGreeness <- function(file_name, clip_coords){ - band_image = raster(file_name, band = 1) - red_crop = crop(band_image, clip_coords) + + band_image <- raster(file_name, band = 1) + red_crop <- crop(band_image, clip_coords) - band_image = raster(file_name, band = 2) - green_crop = crop(band_image, clip_coords) + band_image <- raster(file_name, band = 2) + green_crop <- crop(band_image, clip_coords) add_rasters <- green_crop + red_crop numerator <- cellStats(add_rasters, stat = "sum") @@ -213,12 +214,7 @@ getGreeness <- function(file_name, clip_coords){ denominator <- cellStats(subtract_rasters, stat = "sum") greeness <- numerator / denominator - - # If we can't get a greeness value we return zero - if (is.nan(greeness)) { - greeness = 0; - } - + return(greeness) }