Merge remote-tracking branch 'origin/new_static_downloads' into new_s…

…tatic_downloads

.gitignore

@@ -18,3 +18,8 @@ data/**/*.parquet
 data/**/*.zip
 data/**/*.nc
 data/**/*.tif
+data/elevation/*
+data/slope_aspect/*
+data/soil/*
+
+

R/get_slope_aspect.R

@@ -11,14 +11,17 @@
 #' 
 #' https://www.fao.org/soils-portal/data-hub/soil-maps-and-databases/harmonized-world-soil-database-v12/en/
 library(archive)
-get_slope_aspect<- function(slope_aspect_directory_dataset, slope_aspect_directory_raw, continent_raster_template) {
+
+# NCL: Split into aspect and slope specific functions.
+get_slope_aspect <- function(slope_aspect_directory_dataset, slope_aspect_directory_raw, continent_raster_template) {
 
   slope_aspect <- c("aspect_zero", "aspect_fortyfive", "aspect_onethirtyfive", "aspect_twotwentyfive", "aspect_undef",
                    "slope_zero", "slope_pointfive", "slope_two", "slope_five", "slope_ten", "slope_fifteen", 
                     "slope_thirty", "slope_fortyfive")
 
   for(asp in slope_aspect) { 
 
+    # Change to remove interim file so as not to clog up hd
     url_out<- switch(asp,  "aspect_zero" = "https://www.fao.org/fileadmin/user_upload/soils/HWSD%20Viewer/GloAspectClN_30as.rar",
                     "aspect_fortyfive" = "https://www.fao.org/fileadmin/user_upload/soils/HWSD%20Viewer/GloAspectClE_30as.rar", 
                     "aspect_onethirtyfive" = "https://www.fao.org/fileadmin/user_upload/soils/HWSD%20Viewer/GloAspectClS_30as.rar",
@@ -35,32 +38,38 @@ get_slope_aspect<- function(slope_aspect_directory_dataset, slope_aspect_directo
 
     filename <- paste("data/slope_aspect/", asp, sep="", ".rar")
 
-    #download.file(url=url_out, destfile = filename)
+    download.file(url=url_out, destfile = filename)
 
     rar_name <- file.path(dirname(filename), system2("unrar", c("lb", filename), stdout = TRUE))
     system2("unrar", c("e", "-o+", filename, dirname(filename), ">/dev/null"))
 
+    #
     GloAspectClN_30as <- rast("GloAspectClN_30as.asc")
     GloAspectClE_30as <- rast("GloAspectClE_30as.asc")
     GloAspectClS_30as <- rast("GloAspectClS_30as.asc")
     GloAspectClW_30as <- rast("GloAspectClW_30as.asc")
     GloAspectClU_30as <- rast("GloAspectClU_30as.asc")
 
     raster_stack_aspect <- c(GloAspectClN_30as, GloAspectClE_30as, GloAspectClS_30as, GloAspectClW_30as, GloAspectClU_30as)
+    raster_stack_aspect <- c(GloAspectClN_30as, GloAspectClN_30as + 0.1)
 
-    raster_max_aspect <- max(raster_stack_aspect, na.rm = TRUE)
+    # Which max type logic. Which layer (GloAspectClN_30as, GloAspectClE_30as) has the highest value
+    # Pull out dominant aspect.
+    raster_max_aspect <- which.max(raster_stack_aspect)
 
+    # Transform such that most common aspect is retained
+    # First identify highest percentage aspect per pixel
+    # Then aggregate to most common dominant pixel aspect per cell
     transformed_raster <- transform_raster(raw_raster = rast(raster_max_aspect),
-                                           template = rast(continent_raster_template)) 
-
-
+                                           template = rast(continent_raster_template),
+                                           method = "mode") # Return most common pixel in resampling
 
     #if_else((str_detect(names(transformed_raster), "Aspect", negate=FALSE))==TRUE, c(transformed_raster), 0)
     #if_else statement for aspect and one for slope
 
     # Convert to dataframe
-    #dat_out<- as.data.frame(transformed_raster, xy = TRUE) |> 
-    #  as_tibble() 
+    dat_out<- as.data.frame(transformed_raster, xy = TRUE) |>
+     as_tibble()
 
     #dat_full <- full_join(dat_out, dat_out, by=c("x", "y"))
 

R/preprocess_soil.R

@@ -77,8 +77,7 @@ preprocess_soil <- function(soil_directory_dataset, soil_directory_raw, continen
   #full_join (sand.d1, clay.d1)
 
   #changes from character to factor for the raster
-  for (i in names(records)[c(2:5,7:13,16:17,19:23)])
-  {
+  for (i in names(records)[c(2:5,7:13,16:17,19:23)]) {
     eval(parse(text=paste0("records$",i," <- as.factor(records$",i,")")))
   }
 
@@ -95,13 +94,30 @@ preprocess_soil <- function(soil_directory_dataset, soil_directory_raw, continen
   dat_out <- as.data.frame(hwsd.zhnj.texture, xy = TRUE) |> 
     as_tibble() 
 
-  dat_out$HWSD2 <- if_else(dat_out$HWSD2=="5", "1",
-                   if_else(dat_out$HWSD2=="7", "2",
-                   if_else(dat_out$HWSD2=="8", "3",
-                   if_else(dat_out$HWSD2=="9", "4",
-                   if_else(dat_out$HWSD2=="10", "5",
-                   if_else(dat_out$HWSD2=="11", "6",
-                   if_else(dat_out$HWSD2=="12", "7","0")))))))
+  # At this point:
+  # 1 - clay (heavy)
+  # 2 - silty clay
+  # 3 - clay
+  # 4 - silty clay loam
+  # 5 - clay loam
+  # 6 - silt
+  # 7 - silt loam
+  # 8 - sandy clay
+  # 9 - loam
+  # 10 - sandy clay loam
+  # 11 - sandy loam
+  # 12 - loamy sand
+  # 13 - sand
+
+  # Re-code factor levels to collapse simplex. 
+  # Figure out where key is for the units are in HWSD2
+  dat_out$HWSD2 <- if_else(dat_out$HWSD2=="5", "1", # clay (heavy) + clay loam
+                   if_else(dat_out$HWSD2=="7", "2", # silty clay + silty loam aka
+                   if_else(dat_out$HWSD2=="8", "3", # clay + sandy clay
+                   if_else(dat_out$HWSD2=="9", "4", # silty clay loam
+                   if_else(dat_out$HWSD2=="10", "5", # clay loam + sandy clay loam BUT SEE RULE 1!!!
+                   if_else(dat_out$HWSD2=="11", "6", # silt sandy + loam
+                   if_else(dat_out$HWSD2=="12", "7","0"))))))) # loamy sand + silt loam
 
 
   #create matrix of map unit ids and the variable of interest - DRAINAGE

R/transform_raster.R

@@ -1,11 +1,11 @@
-transform_raster <- function(raw_raster, template) {
+transform_raster <- function(raw_raster, template, method = "cubicspline") {
 
   if(!identical(crs(raw_raster), crs(template))) {
     raw_raster <- terra::project(raw_raster, template)
   }
   if(!identical(origin(raw_raster), origin(template)) ||
      !identical(res(raw_raster), res(template))) {
-    raw_raster <- terra::resample(raw_raster, template, method = "cubicspline")
+    raw_raster <- terra::resample(raw_raster, template, method = method) #
   } 
 
   return(raw_raster)

_targets.R

@@ -90,6 +90,10 @@ tar_target(elevation_directory_dataset,
 tar_target(elevation_preprocessed, 
            process_elevation(elevation_directory_dataset, elevation_downloaded, elevation_directory_raw, continent_raster_template)),
 
+# Any missing static layers?
+# bioclim
+# forest cover
+#
 
 )
 # Dynamic Data Download -----------------------------------------------------------