Downsample (or extract) HWSD data

downsample_extract_hwsd.R

# 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")