...
 
Commits (9)
......@@ -68,13 +68,8 @@ for lead in {1..9} ; do
wgrib2 $TEMP_GRIB -netcdf $OUTFILE_TEMP
rm $TEMP_GRIB
ncrename -h \
-vlatitude,lat \
-vlongitude,lon \
-dlatitude,lat \
-dlongitude,lon \
"-v${GRIB_VAR},${VAR}" \
$OUTFILE_TEMP
ncrename -h -vlatitude,lat -vlongitude,lon ${OUTFILE_TEMP}
ncrename -h "-v${GRIB_VAR},${VAR}" ${OUTFILE_TEMP}
# Drop the time dimension
ncwa -h -O -a time $OUTFILE_TEMP $OUTFILE_TEMP
......
......@@ -50,11 +50,14 @@ wgrib2 $TEMP_GRB2 -nc_grads -netcdf $TEMP_NC1
# command.
ncrename -h -vlatitude,lat $TEMP_NC1
ncrename -h -vlongitude,lon $TEMP_NC1
ncrename -h -vTMP_2maboveground,tmp2m $TEMP_NC1
ncrename -h -vPRATE_surface,prate $TEMP_NC1
ncrename -h -vTMP_2maboveground,T $TEMP_NC1
ncrename -h -vPRATE_surface,Pr $TEMP_NC1
ncrename -h -dlatitude,lat $TEMP_NC1
ncrename -h -dlongitude,lon $TEMP_NC1
ncatted -h -a standard_name,T,c,c,'surface_temperature' $TEMP_NC1
ncatted -h -a standard_name,Pr,c,c,'precipitation_flux' $TEMP_NC1
# Use the --no_tmp_fl option to NCO and manage the temp files ourselves
# Do this to avoid "permission denied" errors when writing to CIFS shares
# under certain conditions
......@@ -71,8 +74,8 @@ ncatted -h -O \
-a longitude_of_prime_meridian,crs,c,d,0 \
-a semi_major_axis,crs,c,d,6378137 \
-a inverse_flattening,crs,c,d,298.257223563 \
-a grid_mapping,tmp2m,c,c,'crs' \
-a grid_mapping,prate,c,c,'crs' \
-a grid_mapping,T,c,c,'crs' \
-a grid_mapping,Pr,c,c,'crs' \
$TEMP_NC4
# Compress the netCDF
......
......@@ -56,6 +56,9 @@ def download(url, output_dir):
sys.stdout.write('.')
sys.stdout.flush()
outfile.write(chunk)
except:
os.remove(os.path.join(output_dir, fname))
raise
finally:
sys.stdout.write('\n')
......@@ -68,29 +71,40 @@ def main(raw_args):
day = int(args.timestamp[6:8])
hour = int(args.timestamp[8:10])
gribfile = "flxf.01.{TIMESTAMP}.{TARGET}.avrg.grib.grb2".format(TIMESTAMP=args.timestamp,
TARGET=args.target)
hindcast = year <= 2009
start_of_rolling_archive = datetime.datetime.now() - datetime.timedelta(days=8) # should have 7 days but could have more or less
timestamp_datetime = datetime.datetime(year, month, day, hour)
if timestamp_datetime > datetime.datetime.utcnow():
print("Can't download forecast with timestamp in the future.", file=sys.stderr)
sys.exit(1)
if hindcast:
grib_pattern = "flxf{TIMESTAMP}.01.{TARGET}.avrg.grb2"
else:
grib_pattern = "flxf.01.{TIMESTAMP}.{TARGET}.avrg.grib.grb2"
url_patterns = []
gribfile = grib_pattern.format(TIMESTAMP=args.timestamp, TARGET=args.target)
if timestamp_datetime > start_of_rolling_archive:
print("Attempting rolling archive URL")
url_patterns += (
'https://nomads.ncep.noaa.gov/pub/data/nccf/com/cfs/prod/cfs/cfs.{YEAR:04d}{MONTH:02d}{DAY:02d}/{HOUR:02d}/monthly_grib_01/{GRIBFILE}',
)
if hindcast:
url_patterns = [
'https://nomads.ncdc.noaa.gov/modeldata/cmd_mm_9mon/{YEAR:04d}/{YEAR:04d}{MONTH:02d}/{YEAR:04d}{MONTH:02d}{DAY:02d}/{GRIBFILE}'
]
else:
print("Attempting long-term archive URL")
url_patterns += (
start_of_rolling_archive = datetime.datetime.now() - datetime.timedelta(days=8) # should have 7 days but could have more or less
timestamp_datetime = datetime.datetime(year, month, day, hour)
if timestamp_datetime > datetime.datetime.utcnow():
print("Can't download forecast with timestamp in the future.", file=sys.stderr)
sys.exit(1)
url_patterns = [
'ftp://nomads.ncdc.noaa.gov/modeldata/cfsv2_forecast_mm_9mon/{YEAR:04d}/{YEAR:04d}{MONTH:02d}/{YEAR:04d}{MONTH:02d}{DAY:02d}/{TIMESTAMP}/{GRIBFILE}',
'https://nomads.ncdc.noaa.gov/modeldata/cfsv2_forecast_mm_9mon/{YEAR:04d}/{YEAR:04d}{MONTH:02d}/{YEAR:04d}{MONTH:02d}{DAY:02d}/{TIMESTAMP}/{GRIBFILE}',
)
]
rolling_url_pattern = 'https://nomads.ncep.noaa.gov/pub/data/nccf/com/cfs/prod/cfs/cfs.{YEAR:04d}{MONTH:02d}{DAY:02d}/{HOUR:02d}/monthly_grib_01/{GRIBFILE}'
if timestamp_datetime > start_of_rolling_archive:
print("Attempting rolling archive URL first")
url_patterns.insert(0, rolling_url_pattern)
else:
print("Attempting long-term archive URL first")
url_patterns.append(rolling_url_pattern)
for url_pattern in url_patterns:
url = url_pattern.format(YEAR=year,
......
#!/usr/bin/env Rscript
# Copyright (c) 2018 ISciences, LLC.
# Copyright (c) 2018-2019 ISciences, LLC.
# All rights reserved.
#
# WSIM is licensed under the Apache License, Version 2.0 (the "License").
......@@ -28,7 +28,8 @@ Options:
'->usage
suppressMessages({
require(wsim.io)
library(wsim.io)
library(wsim.lsm)
})
logging_init('read_binary_grid')
......@@ -151,9 +152,9 @@ standard_attrs <- list(
list(var='T', key='units', val='degree_Celsius')
),
P=list(
list(var='P', key='long_name', val='Precipitation'),
list(var='P', key='standard_name', val='precipitation_amount'),
list(var='P', key='units', val='mm')
list(var='P', key='long_name', val='Precipitation Rate'),
list(var='P', key='standard_name', val='precipitation_flux'),
list(var='P', key='units', val='kg/m^2/s')
)
)
......@@ -176,9 +177,6 @@ main <- function(raw_args) {
}
}
year <- as.integer(substr(args$yearmon, 1, 4))
month <- as.integer(substr(args$yearmon, 5, 6))
infile <- file(args$input, 'rb')
attrs <- standard_attrs[[args$var]]
......@@ -201,6 +199,11 @@ main <- function(raw_args) {
info('Read data for', args$yearmon)
}
if (args$var == 'P') {
# convert mm to mm/s (kg/m^s/s)
dat <- dat / wsim.lsm::days_in_yyyymm(args$yearmon) / 86400
}
to_write <- list()
to_write[[args$var]] <- dat
......
# Copyright (c) 2018 ISciences, LLC.
# Copyright (c) 2018-2019 ISciences, LLC.
# All rights reserved.
#
# WSIM is licensed under the Apache License, Version 2.0 (the "License").
......@@ -27,7 +27,7 @@ read_fits_from_cdf <- function(files) {
extent <- NULL
for (file in files) {
fit <- read_vars_to_cube(file, attrs_to_read=c('distribution', 'variable'))
fit <- read_vars_to_cube(file, attrs_to_read=c('distribution', 'variable', 'units'))
attr(fit, 'filename') <- file
var <- attr(fit, 'variable')
......
# Copyright (c) 2018 ISciences, LLC.
# Copyright (c) 2018-2019 ISciences, LLC.
# All rights reserved.
#
# WSIM is licensed under the Apache License, Version 2.0 (the "License").
......@@ -14,20 +14,62 @@
#' Read model forcing from a netCDF file
#'
#' @param fname netCDF file containing forcing data
#' @param yearmon year and month associated with the forcing, in YYYYMM format
#' @return \code{wsim.lsm.forcing} object containing model forcing
#'
#' @export
read_forcing_from_cdf <- function(fname) {
read_forcing_from_cdf <- function(fname, yearmon) {
contents <- wsim.io::read_vars_from_cdf(fname)
temp_units <- attr(contents$data$T, 'units')
if (temp_units != 'degree_Celsius') {
contents$data$T <- temp_celsius(contents$data$T, temp_units)
attr(contents$data$T, 'units') <- 'degree_Celsius'
}
precip_units <- attr(contents$data$Pr, 'units')
if (precip_units != 'mm') {
contents$data$Pr <- precip_mm(contents$data$Pr, precip_units, days_in_yyyymm(yearmon), 'day')
attr(contents$data$Pr, 'units') <- 'mm'
}
args <- c(contents["extent"],
contents$data[c("pWetDays", "T", "Pr")])
# TODO use udunits2 package to pick out synonyms, or
# perform automatic conversions?
wsim.io::check_units(contents, 'T', 'degree_Celsius', fname)
wsim.io::check_units(contents, 'Pr', 'mm', fname)
return(do.call(make_forcing, args))
}
#' Convert temperature to degrees Celsius
#'
#' @param x temperature value
#' @param u units of \code{x}
#' @return temperature value in degrees Celsius
temp_celsius <- function(x, u) {
units::drop_units(units::set_units(units::set_units(x, u, mode='standard'), 'degree_C'))
}
#' Convert precipitation value to millimeters
#'
#' @param precip precipitation values
#' @param precip_units units of \code{precip}
#' @param duration an optional duration over which precipitation falls
#' @param duration_units units of \code{duration}
#' @return precipitation amount in millimeters
precip_mm <- function(x, u, duration, duration_units) {
water_density <- units::set_units(1000, 'kg/m^3')
# easiest case: we were given precipitation as [L]
try(return(units::drop_units(units::set_units(units::set_units(x, u, mode='standard'), 'mm'))), silent=TRUE)
if (!is.null(duration)) {
duration <- units::set_units(duration, duration_units, mode='standard')
# or maybe we were given a precipitation rate [L/T]
try(return(units::drop_units(units::set_units(units::set_units(x, u, mode='standard') * duration, 'mm'))), silent=TRUE)
# or maybe we were given a mass-based precipitation rate [M/L^2/T]
try(return(units::drop_units(units::set_units(units::set_units(x, u, mode='standard') / water_density * duration, 'mm'))), silent=TRUE)
}
stop(sprintf('Cannot process precipitation data with units of %s', u))
}
......@@ -58,3 +58,13 @@ test_that('coalesce errors out for bad inputs', {
coalesce(dat, matrix(1:6, nrow=2), 'must be a constant or .* same size')
)
})
test_that('we can convert forcing units to those required by the model', {
expect_equal(5, precip_mm(5, 'mm'))
expect_equal(5, precip_mm(0.5, 'cm'))
expect_equal(5, precip_mm(1.929012e-06, 'mm/s', 30, 'days'), tolerance=1e-6)
expect_equal(5, precip_mm(1.929012e-06, 'kg/m^2/s', 30, 'days'), tolerance=1e-6)
expect_equal(5, temp_celsius(5, 'degree_Celsius'))
expect_equal(5, temp_celsius(278.15, 'K'))
})
#!/usr/bin/env Rscript
# Copyright (c) 2018 ISciences, LLC.
# Copyright (c) 2018-2019 ISciences, LLC.
# All rights reserved.
#
# WSIM is licensed under the Apache License, Version 2.0 (the "License").
......@@ -30,6 +30,20 @@ Options:
--append Append output to existing file
'->usage
check_unit_consistency <- function(var, retro_units, forecast_units) {
if (is.null(retro_units)) {
wsim.io::warnf("Unspecified units for retrospective forecast fit of %s", var)
}
if (is.null(forecast_units)) {
wsim.io::warnf("Unspecified units for forecast of %s", var)
}
if (!is.null(retro_units) && !is.null(forecast_units) && retro_units != forecast_units) {
stop(sprintf("Forecast uses units of %s but retrospective forecast distribution uses units of %s.", forecast_units, retro_units))
}
}
main <- function(raw_args) {
args <- wsim.io::parse_args(usage, raw_args)
......@@ -69,6 +83,14 @@ main <- function(raw_args) {
distribution <- attr(retro_fit, 'distribution')
stopifnot(distribution == attr(obs_fit, 'distribution'))
retro_units <- attr(retro_fit, 'units')
forecast_units <- attr(forecast$data[[var]], 'units')
observed_units <- attr(obs_fit, 'units')
# forecast units need not be the same as observed distribution units, since
# we're just matching quantiles
check_unit_consistency(var, retro_units, forecast_units)
corrected[[var]] = wsim.distributions::forecast_correct(distribution,
forecast$data[[var]],
retro_fit,
......@@ -77,12 +99,12 @@ main <- function(raw_args) {
wsim.io::info('Computed bias-corrected values for', var)
attrs <- c(attrs, list(
list(var=var, key="comment", val=paste0("bias-corrected according to ",
distribution,
" fit data in ",
attr(retro_fit, 'filename'),
" and ",
attr(obs_fit, 'filename')))
list(var=var, key="comment", val=sprintf("bias-corrected according to %s fit data in %s and %s",
distribution,
attr(retro_fit, 'filename'),
attr(obs_fit, 'filename'))),
list(var=var, key="units", val=observed_units),
list(var=var, key="standard_name", val=attr(forecast$data[[var]], 'standard_name'))
))
}
}
......
#!/usr/bin/env Rscript
# Copyright (c) 2018 ISciences, LLC.
# Copyright (c) 2018-2019 ISciences, LLC.
# All rights reserved.
#
# WSIM is licensed under the Apache License, Version 2.0 (the "License").
......@@ -79,7 +79,7 @@ main <- function(raw_args) {
for (i in seq_along(forcings)) {
iter_num <- iter_num + 1
forcing <- wsim.lsm::read_forcing_from_cdf(forcings[i])
forcing <- wsim.lsm::read_forcing_from_cdf(forcings[i], state$yearmon)
wsim.io::info("Running LSM for", state$yearmon, "using", forcings[i], "...")
iter <- wsim.lsm::run(static, state, forcing)
......