generated from ecohealthalliance/container-template
-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
5 changed files
with
264 additions
and
38 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) | ||
|
|
||
|
|
||
|
|
||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) | ||
|
|
||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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") | ||
|
|
||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters