diff --git a/R/get_glw_data.R b/R/get_glw_data.R new file mode 100644 index 0000000..1098827 --- /dev/null +++ b/R/get_glw_data.R @@ -0,0 +1,33 @@ +#' .. content for \description{} (no empty lines) .. +#' +#' .. content for \details{} .. +#' +#' @title +#' @param +#' @return +#' @author Whitney Bagge +#' @export +#' +get_glw_data <- function(glw_directory_raw) { + + options(timeout=200) + + location <- c("url_cattle", "url_sheep", "url_goats") + + for(loc in location) { + + url_out<- switch(loc, "url_cattle" = "https://dataverse.harvard.edu/api/access/datafile/6769710", + "url_sheep" = "https://dataverse.harvard.edu/api/access/datafile/6769629", + "url_goats" = "https://dataverse.harvard.edu/api/access/datafile/6769692") + + filename <- paste("data/glw/", loc, sep="", ".tif") + + download.file(url=url_out, destfile = filename) + + } + + return(glw_directory_raw) + + + +} diff --git a/R/preprocess_glw_data.R b/R/preprocess_glw_data.R new file mode 100644 index 0000000..d2d0f08 --- /dev/null +++ b/R/preprocess_glw_data.R @@ -0,0 +1,38 @@ +#' .. content for \description{} (no empty lines) .. +#' +#' .. content for \details{} .. +#' +#' @title +#' @param bounding_boxes +#' @return +#' @author Whitney Bagge +#' @export +#' +library(raster) +library(terra) +preprocess_glw_data<- function(glw_directory_raw, glw_downloaded, continent_raster_template) { + + transformed_raster_cat <- transform_raster(raw_raster = rast(paste0(glw_downloaded, "/url_cattle.tif")), + template = rast(continent_raster_template)) + transformed_raster_sh <- transform_raster(raw_raster = rast(paste0(glw_downloaded, "/url_sheep.tif")), + template = rast(continent_raster_template)) + transformed_raster_go <- transform_raster(raw_raster = rast(paste0(glw_downloaded, "/url_goats.tif")), + template = rast(continent_raster_template)) + + # Convert to dataframe + dat_out_cat<- as.data.frame(transformed_raster_cat, xy = TRUE) |> + as_tibble() + dat_out_sh<- as.data.frame(transformed_raster_sh, xy = TRUE) |> + as_tibble() + dat_out_go<- as.data.frame(transformed_raster_go, xy = TRUE) |> + as_tibble() + + # Save as parquet + write_parquet(dat_out_cat, "data/glw/glw_cattle", compression = "gzip", compression_level = 5) + write_parquet(dat_out_sh, "data/glw/glw_sheep", compression = "gzip", compression_level = 5) + write_parquet(dat_out_go, "data/glw/glw_goats", compression = "gzip", compression_level = 5) + + + return(glw_directory_raw) + +} diff --git a/R/preprocess_soil.R b/R/preprocess_soil.R new file mode 100644 index 0000000..d36c994 --- /dev/null +++ b/R/preprocess_soil.R @@ -0,0 +1,138 @@ +#' .. content for \description{} (no empty lines) .. +#' +#' .. content for \details{} .. +#' +#' @title +#' @param continent_polygon +#' @return +#' @author Whitney Bagge +#' @export +library(DBI) +library(RSQLite) +preprocess_soil <- function(soil_downloaded, continent_raster_template, soil_directory_raw) { + + #read in the raster file + + #crop the raster to the continent + #hwsd_bounded <- terra::crop(unzipped_soil_raster, terra::unwrap(continent_raster_template)) + + #reproject the raster + #print(paste("UTM zone:", utm.zone <- + # floor(((sf::st_bbox(hwsd_bounded)$xmin + + # sf::st_bbox(hwsd_bounded)$xmax)/2 + 180)/6) + # + 1)) + + #(epsg <- 32600 + utm.zone) + + #hwsd_bounded.utm <- project(hwsd_bounded, paste0("EPSG:", epsg), method = "near") + + #terra::resample(hwsd_bounded.utm, method = "near") + + transformed_raster <- transform_raster(raw_raster = rast(paste0(soil_downloaded, "/HWSD2.bil")), + template = rast(continent_raster_template)) + + #connect to database and extract values + m <- dbDriver("SQLite") + con <- dbConnect(m, dbname="data/soil/soil_database.sqlite") + dbListTables(con) + + ####extract map unit codes in bounded area (WINDOW_ZHNJ) to join with SQL databases### + dbWriteTable(con, name="WINDOW_ZHNJ", + value=data.frame(hwsd2_smu = sort(unique(values(transformed_raster)))), + overwrite=TRUE) + + dbExecute(con, "drop table if exists ZHNJ_SMU") # to overwrite + + dbListTables(con) + + #creates a temp database that combines the map unit codes in the raster window to the desired variable + dbExecute(con, + "create TABLE ZHNJ_SMU AS select T.* from HWSD2_SMU as T + join WINDOW_ZHNJ as U + on T.HWSD2_SMU_ID=U.HWSD2_SMU + order by HWSD2_SMU_ID") + + #creates a dataframe "records" in R from SQL temp table created above + records <- dbGetQuery(con, "select * from ZHNJ_SMU") + + #create sand and clay tables in R + #sand.d1 <- dbGetQuery(con, + # "select U.HWSD2_SMU_ID, U.SAND from ZHNJ_SMU as T + # join HWSD2_LAYERS as U on T.HWSD2_SMU_ID=U.HWSD2_SMU_ID + # where U.LAYER='D1' + # order by U.HWSD2_SMU_ID") + # + #clay.d1 <- dbGetQuery(con, + # "select U.HWSD2_SMU_ID, U.CLAY from ZHNJ_SMU as T + # join HWSD2_LAYERS as U on T.HWSD2_SMU_ID=U.HWSD2_SMU_ID + # where U.LAYER='D1' + # order by U.HWSD2_SMU_ID") + + #remove the temp tables and database connection + dbRemoveTable(con, "WINDOW_ZHNJ") + dbRemoveTable(con, "ZHNJ_SMU") + dbDisconnect(con) + + #join sand and clay data frames in r to create a ratio variable + #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)]) + { + eval(parse(text=paste0("records$",i," <- as.factor(records$",i,")"))) + } + + #create matrix of map unit ids and the variable of interest - TEXTURE CLASS + rcl.matrix.texture <- cbind(id = as.numeric(as.character(records$HWSD2_SMU_ID)), + texture = as.numeric(records$TEXTURE_USDA)) + + #classify the raster (transformed_raster) using the matrix of values - TEXTURE CLASS + hwsd.zhnj.texture <- classify(transformed_raster, rcl.matrix.texture) + hwsd.zhnj.texture <- as.factor(hwsd.zhnj.texture) + levels(hwsd.zhnj.texture) <- levels(records$TEXTURE_USDA) + + # Convert to dataframe + 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"))))))) + + + #create matrix of map unit ids and the variable of interest - DRAINAGE + rcl.matrix.drainage <- cbind(id = as.numeric(as.character(records$HWSD2_SMU_ID)), + drainage = as.numeric(records$DRAINAGE)) + + #classify the raster (transformed_raster) using the matrix of values - DRAINAGE + hwsd.zhnj.drainage <- classify(transformed_raster, rcl.matrix.drainage) + hwsd.zhnj.drainage <- as.factor(hwsd.zhnj.drainage) + levels(hwsd.zhnj.drainage) <- levels(records$DRAINAGE) + + # Convert to dataframe + dat_out2 <- as.data.frame(hwsd.zhnj.drainage, xy = TRUE) |> + as_tibble() + + dat_out2$HWSD2 <- if_else(dat_out2$HWSD2=="MW", "4", + if_else(dat_out2$HWSD2=="P", "6", + if_else(dat_out2$HWSD2=="SE", "2", + if_else(dat_out2$HWSD2=="VP", "7","0")))) + + dat_out2$HWSD2 <- as.numeric(as.character(dat_out2$HWSD2)) + + # Save as parquet + write_parquet(dat_out, "data/soil/soil_texture", compression = "gzip", compression_level = 5) + write_parquet(dat_out2, "data/soil/soil_drainage", compression = "gzip", compression_level = 5) + + #writeRaster(hwsd.zhnj.drainage, "data/soil/drainage_raster.tif", overwrite=TRUE) + #writeRaster(hwsd.zhnj.texture, "data/soil/texture_class_raster.tif", overwrite=TRUE) + #writeRaster(x, sand_clay_raster, overwrite=TRUE) + + + return("/data/soil") + +} diff --git a/R/soil_download.R b/R/soil_download.R index 44bcb05..38e577b 100644 --- a/R/soil_download.R +++ b/R/soil_download.R @@ -7,29 +7,30 @@ #' @return #' @author Whitney Bagge #' @export -soil_download <- function() { +soil_download <- function(soil_directory_raw) { options(timeout=200) - loc <- c("soil_database", "soil_raster") - for(loc in location) { + location <- c("soil_database", "soil_raster") + + for(loc in location){ - url_out<- switch(loc, "soil_raster" = "https://s3.eu-west-1.amazonaws.com/data.gaezdev.aws.fao.org/HWSD/HWSD2_RASTER.zip", - "soil_database" = "https://www.isric.org/sites/default/files/HWSD2.sqlite") + url_out<- switch(loc, "soil_raster" = "https://s3.eu-west-1.amazonaws.com/data.gaezdev.aws.fao.org/HWSD/HWSD2_RASTER.zip", + "soil_database" = "https://www.isric.org/sites/default/files/HWSD2.sqlite") - file_ext<- switch(loc,"soil_raster" = ".zip", "soil_database" = ".sqlite") - - file_name <- paste("data/soil/", loc, file_ext, sep="") - - download.file(url=url_out, destfile = file_name) - - if (loc == "soil_raster" ){ - unzip(file_name, exdir = "data/soil/") - } + file_ext<- switch(loc,"soil_raster" = ".zip", "soil_database" = ".sqlite") + + filename <- paste("data/soil/", loc, file_ext, sep="") + + download.file(url=url_out, destfile = filename) + + if (loc == "soil_raster" ){ + unzip(filename, exdir = "data/soil/") + } } - return(file.path(download_directory, filename)) + return(soil_directory_raw) } diff --git a/_targets.R b/_targets.R index f5dd248..b6aae0b 100644 --- a/_targets.R +++ b/_targets.R @@ -38,12 +38,47 @@ static_targets <- tar_plan( tar_target(continent_bounding_box, sf::st_bbox(continent_polygon)), tar_target(continent_raster_template, wrap(terra::rast(ext(continent_polygon), - resolution = 0.5))), #TODO change to 0.1 (might cause error in transform, leaving at 0.5 for now) + resolution = 0.1))), #TODO change to 0.1 (might cause error in transform, leaving at 0.5 for now) # nasa power resolution = 0.5; enmwf = ; ndvi = # tar_target(continent_raster_template_plot, create_raster_template_plot(rast(continent_raster_template), continent_polygon)) - -) + # SOIL ----------------------------------------------------------- + tar_target(soil_directory_raw, + create_data_directory(directory_path = "data/soil")), + tar_target(soil_downloaded, soil_download(soil_directory_raw), + format = "file", + repository = "local"), + tar_target(soil_preprocessed, + preprocess_soil(soil_downloaded, continent_raster_template)), + + # SLOPE and ASPECT ------------------------------------------------- + tar_target(slope_aspect_directory_raw, + create_data_directory(directory_path = "data/slope_aspect")), + tar_target(slope_aspect_downloaded, get_slope_aspect(slope_aspect_directory_raw, continent_polygon), + format = "file", + repository = "local"), + + # Gridded Livestock of the world ----------------------------------------------------------- + tar_target(glw_directory_raw, + create_data_directory(directory_path = "data/glw")), + tar_target(glw_downloaded, get_glw_data(glw_directory_raw), + format = "file", + repository = "local"), + tar_target(glw_preprocessed, + preprocess_glw_data(glw_directory_raw, glw_downloaded, continent_raster_template)), + + +# ELEVATION ----------------------------------------------------------- +tar_target(elevation_directory_raw, + create_data_directory(directory_path = "data/elevation")), +tar_target(elevation_downloaded, get_elevation(elevation_directory_raw), + format = "file", + repository = "local"), +tar_target(elevation_preprocessed, + process_elevation(elevation_directory_raw, elevation_downloaded, continent_raster_template)), + + +) # Dynamic Data Download ----------------------------------------------------------- dynamic_targets <- tar_plan( @@ -52,6 +87,7 @@ dynamic_targets <- tar_plan( tar_target(wahis_rvf_outbreaks_preprocessed, preprocess_wahis_rvf_outbreaks(wahis_rvf_outbreaks_raw)), + # SENTINEL NDVI ----------------------------------------------------------- # 2018-present @@ -239,27 +275,7 @@ dynamic_targets <- tar_plan( format = "file" ), - ## GLW - - tar_target(glw_cattle, get_glw()), - tar_target(glw_sheep, get_glw()), - tar_target(glw_goats, get_glw()), - tar_target(glw_cattle_preprocessed, preprocess_glw(glw_cattle, bounding_boxes)), - #need to cache - - ## ELEVATION - - tar_target(elevation_layer_raw, get_elevation()), - tar_target(elevation_layer_processed, - process_elevation(elevation_layer_processed, bounding_boxes)), - #need to cache - - ## GRAZING CAPACITY - tar_target(grazingcap_raw, get_grazingcap()), - tar_target(grazingcap_processed, - process_grazingcap(grazingcap_layer_processed, bounding_boxes)), - #need to cache )