Commit 2207956e authored by Dan Baston's avatar Dan Baston
Browse files

Merge branch 'r-lsm' into 'master'

Implement LSM in R/C++

See merge request !1
parents f8f72afb 92f1071d
Pipeline #11544654 passed with stage
in 4 minutes and 34 seconds
*.swp
*.pyc
*.o
*.so
RcppExports.cpp
RcppExports.R
.Rproj.user
.idea/
..Rcheck
.Rhistory
......@@ -2,4 +2,4 @@ image: isciences/wsim-gitlabci
check:
script:
- make check
- make install check
install:
for s in wsim.io wsim.distributions wsim.lsm; do \
$(MAKE) -C $${s} $@ || exit 1; \
done;
check:
for s in wsim.io wsim.distributions; do \
for s in wsim.io wsim.distributions wsim.lsm; do \
$(MAKE) -C $${s} $@ || exit 1; \
done;
html:
for s in wsim.io wsim.distributions docs; do \
for s in wsim.io wsim.distributions wsim.lsm docs; do \
$(MAKE) -C $${s} $@ || exit 1; \
done;
......
......@@ -36,3 +36,4 @@ RUN Rscript -e 'install.packages(c("devtools","testthat","roxygen2","raster","Rc
RUN Rscript -e 'devtools::install_github("hadley/pkgdown")'
RUN Rscript -e 'install.packages("foreach", dependencies = TRUE, repos="http://cran.us.r-project.org")'
RUN Rscript -e 'install.packages("lubridate", dependencies = TRUE, repos="http://cran.us.r-project.org")'
wsim_lsm
********
The ``wsim.lsm`` utility runs one ore more iterations of the WSIM Land Surface Model.
Usage is as follows:
.. code-block:: console
wsim_lsm --state <file> \
[--forcing <file>]... \
--flowdir <file> \
--wc <file> \
--elevation <file> \
--results <file> \
--next_state <file>
* ``--state`` is a netCDF file representing the input state of the model. It must provide four variables:
* ``Dr`` amount of detained runoff in millimeters
* ``Ds`` amount of detained snowmelt in millimeters
* ``Snowpack`` snowpack water equivalent in millimeters
* ``Ws`` soil moisture in millimeters
* ``snowmelt_month``` the number of consecutive months of melting conditions
In addition, the state file must define a global attribute ``yearmon`` that specifies a year and month YYYYMM format. The state file will be considered to represent conditions at the start of this month.
* ``--forcing`` specifies one or more netCDF files of forcing data to be used, with each file providing data for a single model iteration. Multiple ``--forcing`` arguments may be provided, and each argument may refer to a single file or a glob of multiple files. Forcing data will be applied in a character-sort order based on the file names of the inputs, not the order in which they are specified. Each forcing file must contain the following variables:
* ``T`` the average monthly temperature in degrees Celsius
* ``Pr`` the total monthly precipitation (rainfall and snowfall) in millimeters
* ``pWetDays`` the fraction of days during which precipitation falls
* ``daylength`` the average day length, as a fraction of a 24-hour period
* ``--wc`` a file providing the soil moisture holding capacity in millimeters
* ``--elevation`` a file providing land surface elevation in meters
* ``--flowdir`` a file providing a surface flow direction matrix (TODO link to specification of format)
The following arguments define model outputs:
* ``--results`` a netCDF file to which model results will be written. If the filename contains the the pattern ``%T``, results from all model iteration will be written to disk as separate files, with the filename formed by substituting the timestep year-month for ``%T``. If the ``%T`` pattern is not present in the filename, only the results of the final iteration will be written to disk.
* ``--next_state`` a netCDF file to which a model state will be written, suitable for use in a subsequent model iteration. Substitution of ``%T`` is performed in the same manner as for the ``--results`` argument.
install:
Rscript -e "devtools::install()"
NAMESPACE:
echo "# Generated by roxygen2: do not edit by hand" > NAMESPACE
check:
check: NAMESPACE
Rscript -e "results <- devtools::check(); stopifnot(length(results[['errors']]) == 0)"
html:
install: NAMESPACE
Rscript -e "devtools::document(); devtools::install()"
html: NAMESPACE
Rscript -e "devtools::document(); pkgdown::build_site()"
check:
NAMESPACE:
echo "# Generated by roxygen2: do not edit by hand" > NAMESPACE
check: NAMESPACE
Rscript -e "results <- devtools::check(); stopifnot(length(results[['errors']]) == 0)"
install:
Rscript -e "devtools::install()"
install: NAMESPACE
Rscript -e "devtools::document(); devtools::install()"
html:
html: NAMESPACE
Rscript -e "devtools::document(); pkgdown::build_site()"
^.*\.Rproj$
^\.Rproj\.user$
Makefile
docs
.Rproj.user
.Rhistory
.RData
docs/
man/
NAMESPACE
Package: wsim.lsm
Title: WSIM Land Surface Model
Version: 0.0.0.9000
Authors@R: person("ISciences, L.L.C.", email = "dbaston@isciences.com", role = c("aut", "cre", "cph"))
Description: The wsim.lsm provides the land surface model used by WSIM.
Depends: R (>= 3.4.0)
License: file LICENSE
Imports:
lubridate,
Rcpp (>= 0.12.12),
raster,
rgdal,
wsim.io
LinkingTo: Rcpp
Suggests:
testthat
Encoding: UTF-8
LazyData: true
RoxygenNote: 6.0.1
NAMESPACE:
echo "# Generated by roxygen2: do not edit by hand" > NAMESPACE
check: NAMESPACE
Rscript -e "results <- devtools::check(); stopifnot(length(results[['errors']]) == 0)"
install: NAMESPACE
Rscript -e "devtools::document(); devtools::install()"
html: NAMESPACE
Rscript -e "devtools::document(); pkgdown::build_site()"
#' @export
cdf_attrs <- list(
list(
var="lat",
description="Latitude (degrees)",
long_name="latitude",
standard_name="latitude",
units="degrees_north"
),
list(
var="lon",
description="Longitude (degrees)",
long_name="longitude",
standard_name="longitude",
units="degrees_north"
),
list(
var="Bt_RO",
description="Total Blue Water (mm)",
long_name="Total Blue Water",
units="mm"
),
list(
var="Bt_Runoff",
description="Total Blue Water, not accounting for detention (mm)",
long_name="Total Blue Water",
units="mm"
),
list(
var="RO_mm",
description="Runoff (mm)",
long_name="Runoff",
standard_name="surface_runoff_amount",
units="mm"
),
list(
var="RO_m3",
description="Runoff (m3)",
long_name="Runoff",
standard_name="surface_runoff_amount",
units="m^3"
),
list(
var="Runoff_mm",
description="Runoff (mm), not accounting for detention",
long_name="Runoff",
standard_name="surface_runoff_amount",
units="mm"
),
list(
var="Runoff_m3",
description="Runoff (m3), not accounting for detention",
long_name="Runoff",
standard_name="surface_runoff_amount",
units="m^3"
),
list(
var="daylength",
description="Length of Day (fraction of 24h period)",
long_name="Daylength",
units="1" # Recommended by CF conventions 1.7, sec 3.1
),
list(
var="pWetDays",
description="Fraction of Days with Precipitation",
long_name="Fraction of Days with Precipitation",
units="1"
),
list(
var="E",
description="Evapotranspiration (mm)",
long_name="Evapotranspiration",
standard_name="water_evaporation_amount",
units="mm"
),
list(
var="EmPET",
description="Evapotranspiration minus Potential Evapotranspiration (mm)",
long_name="Evapotranspiration minus Potential Evapotranspiration",
units="mm"
),
list(
var="PETmE",
description="Potential minus Actual Evapotranspiration (mm)",
long_name="Potential minus Actual Evapotranspiration",
units="mm"
),
list(
var="P_net",
description="Net Precipitation (mm)",
long_name="Net Precipitation",
units="mm"
),
list(
var="Pr",
description="Total Precipitation (mm)",
long_name="Precipitation",
standard_name="precipitation_amount",
units="mm"
),
list(
var="PET",
description="Potential Evapotranspiration (mm)",
long_name="Potential Evapotranspiration",
standard_name="water_potential_evaporation_amount",
units="mm"
),
list(
var="Sm",
description="Snowmelt (mm)",
long_name="Snowmelt",
standard_name="surface_snow_melt_amount",
units="mm"
),
list(
var="Sa",
description="Snow Accumulation (mm)",
long_name="Snow Accmulation",
standard_name="snowfall_amount",
units="mm"
),
list(
var="T",
description="Temperature (C)",
long_name="Temperature",
standard_name="surface_temperature",
units="mm"
),
list(
var="Dr",
description="Detained Runoff (mm)",
long_name="Detained Runoff",
units="mm"
),
list(
var="Ds",
description="Detained Snowmelt (mm)",
long_name="Detained Snowmelt",
units="mm"
),
list(
var="Ws",
description="Soil Moisture (mm)",
long_name="Soil Moisture (mm)",
standard_name="soil_moisture_content",
units="mm"
),
list(
var="Ws_ave",
description="Average Soil Moisture (mm)",
long_name="Average Soil Moisture",
standard_name="soil_moisture_content",
units="mm"
),
list(
var="dWdt",
description="Change in Soil Moisture (mm)",
long_name="Change in Soil Moisture",
units="mm"
),
list(
var="snowmelt_month",
description="Number of consecutive months of melting conditions",
long_name="Snowmelt Month",
units="1"
),
list(
var="Snowpack",
description="Water Equivalent of Snowpack (mm)",
long_name="Snowpack",
standard_name="surface_snow_amount",
units="mm"
)
)
#' Compute the area occupied by each cell of a RasterLayer
#'
#' @param An object that is coercible to a RasterLayer,
#' with a defined extent assumed to be in lon-lat
#' coordinates
#'
#' @return A RasterLayer with the same extent and resolution
#' as the input
#'
#' @export
cell_areas_m2 <- function(rast) {
areas <- raster::raster(rast)
dlon <- raster::res(rast)[1]
dlat <- raster::res(rast)[2]
radius_m <- 6378000
cmp <- sapply(raster::yFromRow(rast, 1:nrow(rast)), function(lat) {
lat1 <- lat - dlat/2
lat2 <- lat + dlat/2
return(pi / 180 * radius_m^2 * abs(sin(lat1 * pi / 180) - sin(lat2 * pi / 180)) * dlon)
})
raster::values(areas) <- rep(cmp, each=ncol(areas))
return(areas)
}
#' Return the number of days in a given month
#'
#' @param yyyymmm Year/month in YYYYMM format
#' @return number of days in month
#' @export
days_in_yyyymm <- function(yyyymm) {
first_day <- as.Date(paste0(yyyymm, '01'), '%Y%m%d')
return(unname(lubridate::days_in_month(first_day)))
}
#' Potential evapotranspiration using Hamon's equation
#'
#' @param lambda Average day length, as a fraction (0 to 1)
#' @param T Average temp (C)
#' @param nDays number of days to estimate
#' @return Potential evapotranspiration (mm)
#' @export
e_potential <- function(lambda, T, nDays) {
# TODO FIXME
# The "237.2" number is a typo, copied from the Kepler workspace
# so that we can verify consistent output.
# The correct number is 273.2
nDays * 715.5 * lambda * saturated_vapor_pressure(T) / (T + 237.2)
}
#' Make a WSIM LSM forcing
#' @export
make_forcing <- function(extent, daylength, pWetDays, T, Pr) {
forcing <- list(
daylength= daylength,
pWetDays= pWetDays,
T= T,
Pr= Pr
)
if (!all(sapply(forcing, is.matrix)))
stop('Non-matrix input in make_forcing')
if (length(unique(lapply(forcing, dim))) > 1)
stop('Unequal matrix dimensions in make_forcing')
forcing$extent <- extent
class(forcing) <- 'wsim.lsm.forcing'
return(forcing)
}
#' Determine if an object represents an LSM forcing
#' @export
is.wsim.lsm.forcing <- function(thing) {
inherits(thing, 'wsim.lsm.forcing')
}
#' Create a set of WSIM LSM results
#' @export
make_results <- function(
Bt_RO,
Bt_Runoff,
E,
EmPET,
PET,
PETmE,
P_net,
RO_m3,
RO_mm,
Runoff_m3,
Runoff_mm,
Sa,
Sm,
Ws_ave,
dWdt,
extent
) {
results <- list(
Bt_RO= Bt_RO,
Bt_Runoff= Bt_Runoff,
E= E,
EmPET= EmPET,
PET= PET,
PETmE= PETmE,
P_net= P_net,
RO_m3= RO_m3,
RO_mm= RO_mm,
Runoff_m3= Runoff_m3,
Runoff_mm= Runoff_mm,
Sa= Sa,
Sm= Sm,
Ws_ave= Ws_ave,
dWdt= dWdt
)
if (!all(sapply(results, is.matrix)))
stop('Non-matrix input in make_results')
if (length(unique(lapply(results, dim))) > 1)
stop('Unequal matrix dimensions in make_results')
results$extent <- extent
class(results) <- "wsim.lsm.results"
return(results)
}
#' Determine if an object represents LSM model results
#' @export
is.wsim.lsm.results <- function(thing) {
inherits(thing, 'wsim.lsm.results')
}
#' Create a WSIM LSM state
#' @export
make_state <- function(extent, Snowpack, Dr, Ds, Ws, snowmelt_month, yearmon) {
matrices <- list(
Snowpack= Snowpack,
Dr= Dr,
Ds= Ds,
Ws= Ws,
snowmelt_month= snowmelt_month
)
attrs <- list(
extent= extent,
yearmon= yearmon
)
if (!all(sapply(matrices, is.matrix)))
stop('Non-matrix input in make_state')
if (length(unique(lapply(matrices, dim))) > 1)
stop('Unequal matrix dimensions in make_state')
if (!(is.character(yearmon) && nchar(yearmon)==6))
stop('Invalid year-month in make_state')
state <- c(matrices, attrs)
class(state) <- 'wsim.lsm.state'
return(state)
}
#' Determine if an object represents an LSM state
#' @export
is.wsim.lsm.state <- function(thing) {
inherits(thing, 'wsim.lsm.state')
}
#' Get the next month
#'
#' @param yyyymmm Year/month in YYYYMM format
#' @return Following month in YYYYMM format
#' @export
next_yyyymm <- function(yyyymm) {
first_day_of_current_month <- as.Date(paste0(yyyymm, '01'), '%Y%m%d')
return(strftime(first_day_of_current_month + 31, '%Y%m'))
}
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