Commit afbd93d7 authored by Cynthia Crowley's avatar Cynthia Crowley

Merge branch 'master' into add-results-vars

parents 65f6eb35 692de30b
Pipeline #76566768 passed with stages
in 17 minutes
#!/usr/bin/env bash
# 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,11 +30,19 @@ TEMP_GRB2=/tmp/regrid_halfdeg.$$.grb2
TEMP_NC1=/tmp/`basename $2`.tmp1.$$.nc
TEMP_NC2=/tmp/`basename $2`.tmp2.$$.nc
TEMP_NC3=/tmp/`basename $2`.tmp3.$$.nc
TEMP_NC4=/tmp/`basename $2`.tmp4.$$.nc
function cleanup {
rm -f $TEMP_GRB2
rm -f $TEMP_NC1
rm -f $TEMP_NC2
rm -f $TEMP_NC3
rm -f $TEMP_NC4
}
trap cleanup EXIT
wgrib2 $1 -match "PRATE:surface|TMP:2 m" -new_grid latlon -179.75:720:0.5 -89.75:360:0.5 $TEMP_GRB2
echo "wgrib2 $TEMP_GRB2 -nc_grads -netcdf $TEMP_NC1"
wgrib2 $TEMP_GRB2 -nc_grads -netcdf $TEMP_NC1
rm $TEMP_GRB2
# Rename each variable in separate commands
# Some versions of NCO fail to find variables
......@@ -56,18 +64,16 @@ ncwa --no_tmp_fl -h -a time $TEMP_NC1 $TEMP_NC2
# Drop the time variable
ncks --no_tmp_fl -h -C -O -x -v time $TEMP_NC2 $TEMP_NC3
# Add a CRS variable
ncap2 --no_tmp_fl -h -O -s 'crs=-9999' $TEMP_NC3 $2
ncap2 --no_tmp_fl -h -O -s 'crs=-9999' $TEMP_NC3 $TEMP_NC4
ncatted -h -O \
-a spatial_ref,crs,c,c,'GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_84\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.017453292519943295]]' \
-a spatial_ref,crs,c,c,'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]' \
-a grid_mapping_name,crs,c,c,'latitude_longitude' \
-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' \
$2
rm $TEMP_NC1
rm $TEMP_NC2
rm $TEMP_NC3
$TEMP_NC4
# Compress the netCDF
ncks $TEMP_NC4 -L1 -7 $2
......@@ -138,15 +138,7 @@ def exact_extract(*,
return Step(
targets=output,
dependencies=[gdaldataset2filename(ds) for ds in rasters.values()] + [boundaries],
commands=[
cmd,
['ncks',
'-O', # overwrite
'-L0', # enable level-0 compression
'-4', # convert to netCDF4. exactextract writes NC3 and then R ncdf4 library will refuse to append to it.
'--fix_rec_dmn', 'id', # convert unlimited id dim to fixed dim. Without this R ncdf4 will not append.
output, output]
],
commands=[cmd],
comment=comment
)
......
#' Create netCDF dimensions
#'
#' @inheritParams write_vars_to_cdf
#' @param unlimited_dims character vector of dimension names which should
#' be created as unlimited
#' @return list of \code{ncdim4} objects
#'
make_netcdf_dims <- function(vars, extent, ids, extra_dims, unlimited_dims=NULL) {
c(make_netcdf_base_dims(vars, extent, ids, unlimited_dims),
make_netcdf_extra_dims(extra_dims, unlimited_dims))
}
make_netcdf_base_dims <- function(vars, extent, ids, unlimited_dims=NULL) {
is_spatial <- is.null(ids)
if (is_spatial) {
lats <- lat_seq(extent, dim(vars[[1]]))
lons <- lon_seq(extent, dim(vars[[1]]))
return(list(
lon = ncdf4::ncdim_def("lon",
units="degrees_east",
vals=lons,
longname="Longitude",
create_dimvar=TRUE,
unlim="lon" %in% unlimited_dims),
lat = ncdf4::ncdim_def("lat",
units="degrees_north",
vals=lats,
longname="Latitude",
create_dimvar=TRUE,
unlim="lat" %in% unlimited_dims)
))
} else {
if (mode(ids) == 'character') {
# The R ncdf4 library does not support proper netCDF 4 strings. So we do it the
# old-school way, with fixed-length character arrays. Data written in this
# way seems to be interpreted correctly by software such as QGIS.
return(list(
id = ncdf4::ncdim_def("id",
units="",
vals=1:length(ids),
create_dimvar=FALSE,
unlim="id" %in% unlimited_dims)
))
} else {
# integer ids
return(list(
id = ncdf4::ncdim_def("id",
units="",
vals=ids,
create_dimvar=TRUE,
unlim="id" %in% unlimited_dims)
))
}
}
}
make_netcdf_extra_dims <- function(extra_dims, unlimited_dims=NULL) {
extra_ncdf_dims <- list()
for (dimname in names(extra_dims)) {
vals <- extra_dims[[dimname]]
unlim <- dimname %in% unlimited_dims
if (mode(vals) == 'character') {
new_dim <- ncdf4::ncdim_def(dimname, units='', vals=1:length(vals), create_dimvar=FALSE, unlim=unlim)
} else {
new_dim <- ncdf4::ncdim_def(dimname, units='', vals=vals, create_dimvar=TRUE, unlim=unlim)
}
extra_ncdf_dims[[dimname]] <- new_dim
}
return(extra_ncdf_dims)
}
......@@ -322,9 +322,13 @@ check_var_list <- function(cdf, vars) {
#' @param dim_values a list whose keys represent dimension names and values
#' represent values along those dimensions
find_offset <- function(cdf, real_dims, dim_values) {
if (is.null(dim_values)) {
return(NA_integer_)
}
sapply(real_dims, function(dimname) {
if (dimname %in% names(dim_values)) {
i <- which(ncdf4::ncvar_get(cdf, dimname)==dim_values[[dimname]])
i <- which(ncdf4::ncvar_get(cdf, dimname) == dim_values[[dimname]])
if (length(i) == 0) {
stop(sprintf("Invalid value \"%s\" for dimension \"%s\".", dim_values[[dimname]], dimname))
}
......@@ -339,5 +343,9 @@ find_offset <- function(cdf, real_dims, dim_values) {
#'
#' @inheritParams find_offset
find_count <- function(real_dims, dim_values) {
if (is.null(dim_values)) {
return(NA_integer_)
}
ifelse(real_dims %in% names(dim_values), 1, -1)
}
# Copyright (c) 2019 ISciences, LLC.
# All rights reserved.
#
# WSIM is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#' Generate a list of standard attributes to include in netCDF outputs
#'
#' @param is_new is the file new?
#' @param is_spatial does the file have spatial data?
#' @param existing_history optional character vector with an existing history
#' entry to which we should potential apppend
#' @return list of attributes
standard_netcdf_attrs <- function(is_new, is_spatial, existing_history=NULL) {
history <- paste0(date_string(), ': ', get_command(), '\n')
if (!is.null(existing_history)) {
if (endsWith(existing_history, history)) {
# skip adding history to avoid duplicate entry
history <- existing_history
} else if (nchar(existing_history) > 0 && !endsWith(existing_history, '\n')) {
# previous history entry didn't end with a newline, so we add it here
history <- paste0(existing_history, '\n', history)
} else {
history <- paste0(existing_history, history)
}
}
ret <- list(
list(key="wsim_version", val=wsim_version_string()),
list(key="history", val=history)
)
if (is_new) {
ret <- c(ret, list(
list(key="date_created", val=date_string())
))
}
if (is_spatial) {
ret <- c(ret, list(
# CF conventions require that each measurement can be located
# on the surface of the Earth. So we only stamp our file with
# a "Conventions" attribute when writing spatial data.
list(key="Conventions", val="CF-1.6"),
list(var="lon", key="axis", val="X"),
list(var="lon", key="standard_name", val="longitude"),
list(var="lat", key="axis", val="Y"),
list(var="lat", key="standard_name", val="latitude")
))
}
return(ret)
}
# Calculate a time stamp on package load so that it's the same for
# each write within a loop. This prevents duplicate history
# entries.
time_loaded <- Sys.time()
date_string <- function() {
strftime(time_loaded, '%Y-%m-%dT%H:%M:%S%z')
}
This diff is collapsed.
# 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").
......@@ -238,7 +238,7 @@ test_that('write_vars_to_cdf provides useful errors if extent is not correctly s
)
expect_error(write_vars_to_cdf(data, fname),
"Must provide either extent or xmin,")
"Must provide either extent or")
expect_error(write_vars_to_cdf(data, fname, extent=c(0, 1, 0, 1), xmin=2),
"Both extent and xmin.* provided")
......@@ -446,6 +446,7 @@ test_that('we can write data where different variables have different numbers of
append=TRUE)
file.remove(fname)
succeed()
})
test_that('we can write multidimensional id-based data', {
......@@ -518,3 +519,43 @@ test_that('We can pass an array that has dimnames instead of a named list of var
file.remove(fname)
})
test_that('we can append to a file created with an unlimited id dim', {
fname <- tempfile(fileext='.nc')
ids <- 3:17
precip <- runif(length(ids))*10
runoff <- precip*runif(length(ids))
id_dim <- ncdf4::ncdim_def('id', units='', vals=ids, unlim=TRUE)
precip_var <- ncdf4::ncvar_def('precip', units='mm', dim=id_dim, compression=1)
cdf <- ncdf4::nc_create(fname, precip_var)
ncdf4::ncvar_put(cdf, precip_var, precip)
ncdf4::nc_close(cdf)
write_vars_to_cdf(list(runoff=runoff), fname, ids=ids, append=TRUE)
file.remove(fname)
succeed()
})
test_that('we can append to a netcdf version 3 file', {
fname <- tempfile(fileext='.nc')
ids <- 3:17
precip <- runif(length(ids))*10
runoff <- precip*runif(length(ids))
id_dim <- ncdf4::ncdim_def('id', units='', vals=ids)
precip_var <- ncdf4::ncvar_def('precip', units='mm', dim=id_dim)
cdf <- ncdf4::nc_create(fname, precip_var)
ncdf4::ncvar_put(cdf, precip_var, precip)
ncdf4::nc_close(cdf)
write_vars_to_cdf(list(runoff=runoff), fname, ids=ids, append=TRUE)
file.remove(fname)
succeed()
})
# Copyright (c) 2019 ISciences, LLC.
# All rights reserved.
#
# WSIM is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
require(testthat)
context("netCDF attributes")
get_attribute <- function(attrs, attr_name) {
Filter(function(attr) attr$key == attr_name, attrs)[[1]]$val
}
test_that('timestamp is constant between multiple calls', {
first <- date_string()
Sys.sleep(2)
second <- date_string()
expect_equal(first, second)
})
test_that('history entries are not duplicated', {
first <- get_attribute(standard_netcdf_attrs(is_new=FALSE, is_spatial=FALSE), 'history')
second <- get_attribute(standard_netcdf_attrs(is_new=FALSE, is_spatial=FALSE, existing_history=first), 'history')
expect_equal(first, second)
})
......@@ -91,6 +91,7 @@ main <- function(raw_args) {
args$output,
extent= extent,
attrs=attrs,
prec='single',
append=args$append)
wsim.io::info('Wrote corrected forecast to', args$output)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment