Merge pull request #5 from steeleb/main

Small bug fixes

_targets.R

@@ -74,7 +74,8 @@ list(
     name = locs,
     command = locs_save,
     read = read_csv(!!.x),
-    packages = "readr"
+    packages = "readr",
+    error = "null"
   ),
 
   # use location shapefile and configurations to get polygons from NHDPlusv2
@@ -182,8 +183,7 @@ list(
       make_collated_data_files
       add_metadata(yaml = yml,
                    file_prefix = yml$proj,
-                   version_identifier = yml$run_date,
-                   collation_identifier = "2024-08-01")
+                   version_identifier = yml$run_date)
     },
     packages = c("tidyverse", "feather")
   )

config.yml

@@ -9,7 +9,7 @@ local_settings:
 # information is stored - this path must end with a '/'
 - location_file: "" # name of the *.csv* file that contains the location information
 - unique_id: "" # this is the column that stores the unique identifier for each 
-# site, should not contain any special characters
+# site or polygon and should not contain any special characters
 - latitude: "" # this is the column that stores the latitude of the site, must 
 # be in decimal degrees
 - longitude: "" # this is the column that stores the longitude of the site, must 
@@ -43,7 +43,7 @@ temporal_settings:
 # is used, the date will be set to the current date
 
 spatial_settings: 
-- extent: "site" # options: "site", "polygon", "polycenter", "site+poly", 
+- extent: "site" # options: "site", "polygon", "polycenter", "site+polygon", 
 # "site+polygon+polycenter", "polygon+polycenter" - at this time lake and lake 
 # center can only be calculated for lakes in the US
 - site_buffer: 120 # buffer distance in meters around the site or poly center

data_acquisition/py/runGEEperTile.py

@@ -1,6 +1,7 @@
 #import modules
 import ee
 import time
+import fiona
 from datetime import date, datetime
 import os 
 from pandas import read_csv
@@ -1080,7 +1081,7 @@ def ref_pull_89_DSWE3(image):
   Returns:
       summaries for band data within any given geometry area where the DSWE value is 3
   """
-# process image with the radsat mask
+  # process image with the radsat mask
   r = add_rad_mask(image).select('radsat')
   # process image with cfmask
   f = cf_mask(image).select('cfmask')
@@ -1092,7 +1093,7 @@ def ref_pull_89_DSWE3(image):
   h = calc_hill_shades(image, wrs.geometry()).select('hillShade')
   #calculate hillshadow
   hs = calc_hill_shadows(image, wrs.geometry()).select('hillShadow')
-
+  
   #apply dswe function
   d = DSWE(image).select('dswe')
   gt0 = (d.gt(0).rename('dswe_gt0')
@@ -1265,7 +1266,7 @@ def maximum_no_of_tasks(MaxNActive, waitingPeriod):
   locs_feature = csv_to_eeFeat(locations, yml['location_crs'][0])
 
 
-if 'poly' in extent:
+if 'polygon' in extent:
   #if polygon is in extent, check for shapefile
   shapefile = yml['polygon'][0]
   # if shapefile provided by user 
@@ -1287,13 +1288,19 @@ def maximum_no_of_tasks(MaxNActive, waitingPeriod):
   poly_feat = ee.FeatureCollection(features)
 
 
-if 'center' in extent:
+if 'polycenter' in extent:
   if yml['polygon'][0] == True:
     centers_csv = read_csv('data_acquisition/out/user_polygon_centers.csv')
+    centers_csv = (centers_csv.rename(columns={'poi_latitude': 'Latitude', 
+      'poi_longitude': 'Longitude',
+      'r_id': 'id'}))
     # load the shapefile into a Fiona object
-    centers = csv_to_eeFeat(centers_csv, yml['poly_crs'][0])
+    centers = csv_to_eeFeat(centers_csv, 'EPSG:4326')
   else: #otherwise use the NHDPlus file
     centers_csv = read_csv('data_acquisition/out/NHDPlus_polygon_centers.csv')
+    centers_csv = (centers_csv.rename(columns={'poi_latitude': 'Latitude', 
+      'poi_longitude': 'Longitude',
+      'r_id': 'id'}))
     centers = csv_to_eeFeat(centers_csv, 'EPSG:4326')
   # Create an ee.FeatureCollection from the ee.Features
   ee_centers = ee.FeatureCollection(centers)    
@@ -1388,16 +1395,17 @@ def maximum_no_of_tasks(MaxNActive, waitingPeriod):
     feat = (locs_feature
       .filterBounds(geo)
       .map(dp_buff))
-    if e == 'poly':
-      ## get the polygon stack ##
-      feat = (poly_feat
-        .filterBounds(geo))
-      if e == 'center':
-        ## get centers feature and buffer ##
-        feat = (ee_centers
-          .filterBounds(geo)
-          .map(dp_buff))
-      else: print('Extent not identified. Check configuration file.')
+  elif e == 'polygon':
+    ## get the polygon stack ##
+    feat = (poly_feat
+      .filterBounds(geo))
+  elif e == 'polycenter':
+    ## get centers feature and buffer ##
+    feat = (ee_centers
+      .filterBounds(geo)
+      .map(dp_buff))
+  else: 
+    print('Extent not identified. Check configuration file.')
 
   ## process 457 stack
   #snip the ls data by the geometry of the location points    
@@ -1541,16 +1549,17 @@ def maximum_no_of_tasks(MaxNActive, waitingPeriod):
     feat = (locs_feature
       .filterBounds(geo)
       .map(dp_buff))
-    if e == 'poly':
-      ## get the polygon stack ##
-      feat = (poly_feat
-        .filterBounds(geo))
-      if e == 'center':
-        ## get centers feature and buffer ##
-        feat = (ee_centers
-          .filterBounds(geo)
-          .map(dp_buff))
-      else: print('Extent not identified. Check configuration file.')
+  elif e == 'polygon':
+    ## get the polygon stack ##
+    feat = (poly_feat
+      .filterBounds(geo))
+  elif e == 'polycenter':
+    ## get centers feature and buffer ##
+    feat = (ee_centers
+      .filterBounds(geo)
+      .map(dp_buff))
+  else: 
+    print('Extent not identified. Check configuration file.')
 
   # snip the ls data by the geometry of the location points    
   locs_stack_ls89 = (ls89
@@ -1615,7 +1624,7 @@ def maximum_no_of_tasks(MaxNActive, waitingPeriod):
       locs_dataOut_89_D1a.start()
       print('Completed Landsat 8, 9 DSWE 1a stack acquisitions for ' + e + ' configuration at tile ' + str(tiles))
     else:
-      print("Not configured to acquire DSWE 1a stack for Landsat 8, 9 for ' + e + ' configuration")
+      print('Not configured to acquire DSWE 1a stack for Landsat 8, 9 for ' + e + ' configuration')
       print('Starting Landsat 8, 9 DSWE1 acquisition for ' + e + ' configuration at tile ' + str(tiles))
       locs_out_89_D1 = locs_stack_ls89.map(ref_pull_89_DSWE1).flatten()
       locs_out_89_D1 = locs_out_89_D1.filter(ee.Filter.notNull(['med_Blue']))

data_acquisition/src/add_metadata.R

@@ -8,20 +8,17 @@
 #' @param file_prefix specified string that matches the file group to collate
 #' @param version_identifier user-specified string to identify the RS pull these
 #' data are associated with
-#' @param collation_identifier user-specified string to identify the output of this
-#' target
 #' @returns silently creates collated .feather files from 'mid' folder and 
 #' dumps into 'data'
 #' 
 #' 
 add_metadata <- function(yaml,
                          file_prefix, 
-                         version_identifier,
-                         collation_identifier) {
+                         version_identifier) {
 
   files <- list.files(file.path("data_acquisition/mid/"),
-                     pattern = file_prefix,
-                     full.names = TRUE) %>% 
+                      pattern = file_prefix,
+                      full.names = TRUE) %>% 
     # and grab the right version
     .[grepl(version_identifier, .)]
 
@@ -64,20 +61,31 @@ add_metadata <- function(yaml,
     } else {
       ext <- e
     }
-
+    
     # get file using ext
     file <- files[grepl(ext, files)]
     # load file
-    df <- read_feather(file)
+    df <- read_feather(file) %>% 
+      mutate(mission = case_when(grepl("LT04", `system:index`) ~ "LANDSAT_4",
+                                 grepl("LT05", `system:index`) ~ "LANDSAT_5",
+                                 grepl("LE07", `system:index`) ~ "LANDSAT_7",
+                                 grepl("LC08", `system:index`) ~ "LANDSAT_8",
+                                 grepl("LC09", `system:index`) ~ "LANDSAT_9",
+                                 TRUE ~ NA_character_)) 
 
     if (e == "site") {
       spatial_info <- read_csv(file.path(yaml$data_dir,
                                          yaml$location_file)) %>% 
         rename(r_id = yaml$unique_id)%>% 
         mutate(r_id = as.character(r_id))
     } else if (e == "polycenter") {
-      spatial_info <- read_csv("data_acquisition/out/NHDPlus_polygon_centers.csv") %>% 
-        mutate(r_id = as.character(r_id))
+      if (yaml$polygon) { 
+        spatial_info <- read_csv("data_acquisition/out/user_polygon_withrowid.csv") %>% 
+          mutate(r_id = as.character(r_id))
+      } else {
+        spatial_info <- read_csv("data_acquisition/out/NHDPlus_polygon_centers.csv") %>% 
+          mutate(r_id = as.character(r_id))
+      }
     } else if (e == "polygon") {
       if (yaml$polygon) {
         spatial_info <- read_sf(file.path(yaml$poly_dir,
@@ -89,26 +97,41 @@ add_metadata <- function(yaml,
           mutate(r_id = as.character(r_id))
       }
     }
+
     # format system index for join - right now it has a rowid and the unique LS id
     # could also do this rowwise, but this method is a little faster
     df$r_id <- map_chr(.x = df$`system:index`, 
-                            function(.x) {
-                              parsed <- str_split(.x, '_')
-                              str_len <- length(unlist(parsed))
-                              unlist(parsed)[str_len]
-                            })
+                       function(.x) {
+                         parsed <- str_split(.x, '_')
+                         last(unlist(parsed))
+                       })
     df$system.index <- map_chr(.x = df$`system:index`, 
-                                   #function to grab the system index
-                                   function(.x) {
-                                     parsed <- str_split(.x, '_')
-                                     str_len <- length(unlist(parsed))
-                                     parsed_sub <- unlist(parsed)[1:(str_len-1)]
-                                     str_flatten(parsed_sub, collapse = '_')
-                                     })
+                               #function to grab the system index
+                               function(.x) {
+                                 parsed <- str_split(.x, '_')
+                                 str_len <- length(unlist(parsed))
+                                 parsed_sub <- unlist(parsed)[1:(str_len-1)]
+                                 str_flatten(parsed_sub, collapse = '_')
+                               })
+
+    # dswe info is stored differently in each mission group because of character length
+    # so grab out mission-specific dswe info and use that to define dswe
+    mission_dswe <- df %>% 
+      group_by(mission) %>% 
+      slice(1) %>% 
+      ungroup()
+    dswe_loc <- as_tibble(str_locate(mission_dswe$source, "DSWE")) %>% 
+      rowid_to_column() %>% 
+      left_join(., mission_dswe %>% rowid_to_column()) %>% 
+      select(rowid, mission, start, end) %>% 
+      mutate(end = end + 2)
+
     df <- df %>% 
       select(-`system:index`) %>% 
       left_join(., metadata_light) %>% 
-      mutate(DSWE = str_split(source, "_")[[1]][7], .by = source) %>% 
+      left_join(., dswe_loc) %>% 
+      mutate(DSWE = str_sub(source, start, end), .by = mission,
+             DSWE = str_remove(DSWE, "_")) %>%
       left_join(., spatial_info)
 
     # break out the DSWE 1 data
@@ -121,7 +144,7 @@ add_metadata <- function(yaml,
                                      "_collated_DSWE1_",
                                      ext,
                                      "_meta_v",
-                                     collation_identifier,
+                                     version_identifier,
                                      ".feather")))
     } 
 
@@ -135,7 +158,7 @@ add_metadata <- function(yaml,
                                      "_collated_DSWE1a_",
                                      ext, 
                                      "_meta_v",
-                                     collation_identifier,
+                                     version_identifier,
                                      ".feather")))
     }
 
@@ -149,7 +172,7 @@ add_metadata <- function(yaml,
                                      "_collated_DSWE3_",
                                      ext,
                                      "_meta_v",
-                                     collation_identifier,
+                                     version_identifier,
                                      ".feather")))
     }
   })
@@ -159,7 +182,7 @@ add_metadata <- function(yaml,
              pattern = file_prefix,
              full.names = TRUE) %>% 
     #but make sure they are the specified version
-    .[grepl(collation_identifier, .)] %>% 
+    .[grepl(version_identifier, .)] %>% 
     .[!grepl('filtered', .)]
 
 }

data_acquisition/src/calc_center.R

@@ -35,6 +35,10 @@ calc_center <- function(poly, yaml) {
       # Xiao Yang's code in EE - Yang, Xiao. (2020). Deepest point calculation 
       # for any given polygon using Google Earth Engine JavaScript API 
       # (Version v1). Zenodo. https://doi.org/10.5281/zenodo.4136755
+      # this needs to be in WGS for best results
+      if (yaml$poly_crs != "EPSG:4326" & !is.na(yaml$poly_crs)) {
+        one_wbd <- st_transform(one_wbd, "EPSG:4326")
+      }
       coord_for_UTM <- one_wbd %>% st_coordinates()
       mean_x <- mean(coord_for_UTM[ ,1])
       mean_y <- mean(coord_for_UTM[ ,2])