Reading HRLDAS output in R

[fdesales]

Hi. Any recommendations on how to read HRLDAS output files in R? I am able to read WRF output files using the eixport package, however, HRLDAS output files do not contain latitude/longitude (e.g., XLAT_M/XLONG_M, etc) grid information. Thanks!

[cenlinhe]

@RAbolafiaRosenzweig Can you help with this since you are familiar with R?

[PGhaneei]

Hi,

This is an example of how to read and plot outputs from LDASOUT, as well as the grid information from the HRLDAS_setup file in R:

Install packages

install.packages("ncdf4")
install.packages("reshape2")
install.packages("ggplot2")

library(ncdf4)
library(ggplot2)
library(reshape2)

Read output

file_path <- "./2015040300.LDASOUT_DOMAIN1"
nc <- nc_open(file_path)

Extract the soil moisture

smc <- ncvar_get(nc, "SOIL_M")

nc_close(nc)

Select the first layer of sil moisture:

smc_slice <- smc[ , 1, ]

smc_melt <- melt(smc_slice)

Plot

ggplot(smc_melt, aes(Var1, Var2, fill = value)) +
geom_raster() +
scale_fill_gradientn(colors = terrain.colors(10), limits = c(0, 1)) +
theme_minimal() +
labs(title = "Soil Moisture (SOIL_M)", x = "Longitude Index", y = "Latitude Index", fill = "SOIL_M") +
coord_fixed()

Lat and Lon of grids

file_path_grids <- "./HRLDAS_setup_2015040100_d1"
nc_grids <- nc_open(file_path_grids)

Extract the XLAT

XLAT <- ncvar_get(nc_grids, "XLAT")

Extract the XLONG

XLONG <- ncvar_get(nc_grids, "XLONG")

[fdesales]

Thanks @PGhaneei for your script. I was wondering if your ggplot map is correctly projected? In my case, I can get ggplot to generate a map but its geographical projection is wrong. Would you know how to use XLAT and XLONG to reproject the map?

[fdesales]

For example, this is what I get from ggplot using your script. Lat/longs are incorrect. My domain should be center over Botswana. Having the correct projection is crucial for calculating spatial averages and plotting geo boundaries. Thanks

[fdesales]

I have been able to fix the lat/lon issue using the ext function. Please see below. I still to make sure the projection is correct, though. I am using the R terra package to handle raster data more easily.

#Read latitude and longitude from the setup file (forcing)
grd = rast("HRLDAS_setup_2020090100_d1")
xlat = grd[[1]]
xlon = grd[[2]]

#Make a domain with max and min lats and lons
ymn = min(values(xlat))
ymx = max(values(xlat))
xmn = min(values(xlon))
xmx = max(values(xlon))
domain = ext(xmn,xmx,ymn,ymx)

#Read model output and list variable names
dat = rast("2020090500.LDASOUT_DOMAIN1")
names(dat)

#Extract the first variable (IVGTYP) and fix lat/lon
vtype = dat[[1]]
ext(vtype) <- domain  #correct lat/lon

[PGhaneei]

For spatial analysis, what I do is create a new data frame with XLAT and XLONG. Here is an example:

Install necessary packages

install.packages("ncdf4")
install.packages("reshape2")
install.packages("ggplot2")
install.packages("spatstat")
install.packages("sf")
install.packages("sp")

library(ncdf4)
library(ggplot2)
library(reshape2)
library(spatstat)
library(sf)
library(sp)

Read output

file_path <- "./2015040300.LDASOUT_DOMAIN1"
nc <- nc_open(file_path)

Extract the soil moisture

smc <- ncvar_get(nc, "SOIL_M")

nc_close(nc)

Select the first layer of soil moisture:

smc_slice <- smc[ , 1, ]

############## Lat and Lon of grids ############
file_path_grids <- "./HRLDAS_setup_2015040100_d1"
nc_grids <- nc_open(file_path_grids)

Extract the XLAT

XLAT <- ncvar_get(nc_grids, "XLAT")

Extract the XLONG

XLONG <- ncvar_get(nc_grids, "XLONG")

nc_close(nc_grids)

############## Create Complete Data Frame for Spatial Analysis ############

spatial_data <- data.frame(
lon = as.vector(XLONG),
lat = as.vector(XLAT),
soil_m = as.vector(smc_slice)
)

############## Convert to Geo Data Frames ############
spatial_data_sf <- st_as_sf(spatial_data, coords = c("lon", "lat"), crs = 4326)

Plot

ggplot(data = spatial_data_sf) +
geom_sf(aes(color = soil_m)) +
scale_color_viridis_c( limits = c(0, 1)) +
labs(title = "soil_m",
color = "soil_m") +
theme_minimal()

[FuYingWang]

There is a package named 'rwrfhydro' in R.

[FuYingWang]

For multi-layer variables such as' SOIL_M', 'inLyr=1' is also required to extract the value of a certain layer

[fdesales]

I tried using the rwrfhydro package but it did not work for HRLDAS output files.

[FuYingWang]

Thanks @PGhaneei , it help a lot. Is there a way to output 'spatial_data_stf' in tiff format?