Merge pull request #113 from terraref/develop

Remove gdal install requirement

Author: Chris-Schnaufer

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,61 +186,43 @@ 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){
   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)
-
-  # 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)
 }
 
 # Extract all the dates from the images
 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))
 ```