Last active
July 7, 2021 17:38
-
-
Save khufkens/4f9c49ec4ec0341f7b293331b3a56f7c to your computer and use it in GitHub Desktop.
Downsample (or extract) HWSD data
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
# Downsample HWSD data to the scale of 0.5 degree | |
# WFDEI / CRU data | |
library(raster) | |
library(hwsdr) | |
library(tidyverse) | |
# get mu_global at 0.5 degrees from the | |
# ORNL DAAC (which is easier than reading in | |
# that monstrous BIL file) | |
mu_global <- ws_subset( | |
site = "HWSD", | |
location = c(-90, -180, 90, 180), | |
param = "MU_GLOBAL", | |
path = tempdir(), | |
internal = TRUE | |
) | |
# aggregate to 0.5 degrees (comment out if you want ORNL DAAC resolution) | |
mu_global <- aggregate(mu_global, fact = 10, fun = median) | |
# original bil file can be used to the same effect | |
# but needs custom aggretation | |
#mu_global <- raster("hwsd.bil") | |
# select unique values (profiles) to extract | |
# from the SQL database | |
vals <- mu_global %>% | |
raster::unique() | |
vals <- vals[which(vals != 0)] | |
# create an empty table with mu_global index values | |
# grab the mysql database here: | |
# https://github.com/bluegreen-labs/leafnp/raw/main/data-raw/HWSD.sqlite | |
# or elsewhere on github (and fix the path below) | |
con <- DBI::dbConnect(RSQLite::SQLite(), | |
dbname = "data-raw/HWSD.sqlite") | |
DBI::dbWriteTable(con, | |
name = "WINDOW_TMP", | |
value = data.frame(smu_id = vals), | |
overwrite = TRUE | |
) | |
# the real database links indexed locations to "fixed profiles" | |
# using a SQL call and merging the tables | |
result <- DBI::dbGetQuery( | |
con, "select T.* from HWSD_DATA as T join | |
WINDOW_TMP as U on T.mu_global=u.smu_id order by su_sym90") | |
DBI::dbRemoveTable(con, "WINDOW_TMP") | |
# select dominant profile (1) | |
result <- result %>% | |
filter(SEQ == 1) | |
# list all layers in the product | |
layers <- names(result) | |
# loop over all layers and substitute the | |
v <- lapply(layers, function(x){ | |
# subsitute values matching index (id) with v (various result layers) | |
subs(mu_global, data.frame(id = result$MU_GLOBAL, v=result[x])) | |
}) | |
# stack stuff | |
v <- do.call("brick",v) | |
# write values to disk as geotiff | |
writeRaster(v, "data/hwsd_0.5deg.tif",options = c("COMPRESS=DEFLATE")) | |
# convert to dataframe | |
df <- v %>% | |
rasterToPoints() %>% | |
as.data.frame() %>% | |
rename( | |
'latitude' = 'y', | |
'longitude' = 'x' | |
) | |
# write values to disk as rds dataframe | |
saveRDS(df, "data/hwsd_0.5deg.rds") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment