Merge branch 'urbanVPRM'

DESCRIPTION

@@ -1,19 +1,19 @@
 Package: VPRMLandSfcModel
 Type: Package
 Title: R implementation of VPRM with parameter estimation
-Version: 1.2.1
-Date: 2017-10-03
+Version: 1.3.0
+Date: 2023-07-07
 Author: Timothy W. Hilton
-Maintainer: Timothy W. Hilton <[email protected]>
+Maintainer: Timothy W. Hilton <[email protected]>
 Description: Provides an R implementation of the Vegetation
-    Photosynthesis and Respiration Model (VPRM) and tools to estimate
-    VPRM parameter values from observations.  VPRM is a simple
-    diagnostic land surface model based on light-use efficiency.  VPRM
-    diagnoses gross ecosystem exchange (GEE) of carbon dioxide,
-    ecosystem respiration (R), and net ecosystem exchange (NEE) of
-    carbon dioxide. To cite the VPRMLandSfcModel R package or the VPRM
-    model itself in a publication, or for a more detailed description
-    of VPRM, see (within R) 'citation("VPRMLandSurfaceModel")'.
+    Photosynthesis and Respiration Model (VPRM), the urbanVPRM, and tools to
+    estimate VPRM parameter values from observations. VPRM is a simple
+    diagnostic land surface model based on light-use efficiency. VPRM diagnoses
+    gross ecosystem exchange (GEE) of carbon dioxide, ecosystem respiration (R),
+    and net ecosystem exchange (NEE) of carbon dioxide. To cite the
+    VPRMLandSfcModel R package or the (urban)VPRM model itself in a publication, or for
+    a more detailed description of VPRM, see (within R)
+    'citation("VPRMLandSurfaceModel")'.
 License: GPL-3
 LazyLoad: yes
 Imports:
@@ -24,3 +24,5 @@ Suggests:
     knitr,
     rmarkdown
 VignetteBuilder: knitr
+RoxygenNote: 7.2.3
+Encoding: UTF-8

NAMESPACE

@@ -1,11 +1,13 @@
-# Generated by roxygen2 (4.0.2): do not edit by hand
+# Generated by roxygen2: do not edit by hand
 
 S3method(as.data.frame,VPRM_driver_data)
+export(SWDN_2_PAR)
 export(VPRM_driver_data)
 export(detect_large_greenness_change_periods)
 export(estimate_VPRM_pars)
 export(getLSWI)
 export(getPscale)
+export(getPscale_urban)
 export(getTscale)
 export(getWscale)
 export(interpMODIS)

NEWS

@@ -1,3 +1,6 @@
+New features in version 1.3.0
+- implement the urbanVPRM adaptation of VPRM to urban settings
+
 New features in version 1.2.1
 - added vignette to provide usage examples
 
@@ -15,4 +18,4 @@ New features in version 1.1
 - added Park_Falls dataset
 - added example code for all functions
 - changed package name from VPRMModel to VPRMLandSfcModel.  See
-  ChangeLog for reasoning.
+  ChangeLog for reasoning.

R/VPRMLandSfcModel-package.R

@@ -0,0 +1,6 @@
+#' @keywords internal
+"_PACKAGE"
+
+## usethis namespace: start
+## usethis namespace: end
+NULL

R/VPRM_driver_data.R

@@ -5,61 +5,70 @@
 ##' necessary to run VPRM for a single eddy covariance site,
 ##' calculates "derived" fields, and interpolates phenology dynamics.
 ##'
-##' "Derived fields" denote fields that are calculated from other
-##' observed quanitities.  For example, land surface water index
-##' (LSWI) is a derived field, as it is calculated from MODIS
-##' reflectances in the short infrared and near infrared bands.
+##' "Derived fields" denote fields that are calculated from other observed
+##' quanitities. For example, land surface water index (LSWI) is a derived
+##' field, as it is calculated from MODIS reflectances in the short infrared and
+##' near infrared bands.
 ##' @title build a VPRM_driver_data object
-##' @param name_long character string; "short name" of the site.  e.g. US-PFa
-##' @param name_short character string; "long name" of the site.  e.g. Park Falls
+##' @param name_long character string; "short name" of the site. e.g. US-PFa
+##' @param name_short character string; "long name" of the site. e.g. Park Falls
 ##' @param lat numeric; latitude of the site (deg N)
-##' @param lon  numeric; longitude of the site (deg E)
-##' @param PFT character string; plant functional type.  Will be converted to a factor.
-##' @param note character string; optional note; could be anything the
-##' user finds useful.
-##' @param Tmin numeric; minimum temperature for photosynthesis (deg
-##' C).  See Mahadevan et al (2008) eq. 6.
-##' @param Tmax  numeric; maximim temperature for photosynthesis (deg
-##' C).  See Mahadevan et al (2008) eq. 6.
-##' @param Topt  numeric; optimum temperature for photosynthesis (deg
-##' C).  See Mahadevan et al (2008) eq. 6.
-##' @param Tlow numeric; minimum temperature for respiration (deg C).
-##' See Mahadevan et al (2008) eq. 10.
-##' @param tower_date chron vector; timestamps for all tower
-##' observations (NEE_obs, T, PAR)
-##' @param NEE_obs numeric vector; eddy covariance observed net
-##' ecosystem exchange (NEE, umol m-2 s-1)
+##' @param lon numeric; longitude of the site (deg E)
+##' @param PFT character string; plant functional type. Will be converted to a
+##'   factor.
+##' @param note character string; optional note; could be anything the user
+##'   finds useful.
+##' @param Tmin numeric; minimum temperature for photosynthesis (deg C). See
+##'   Mahadevan et al (2008) eq. 6.
+##' @param Tmax numeric; maximim temperature for photosynthesis (deg C). See
+##'   Mahadevan et al (2008) eq. 6.
+##' @param Topt numeric; optimum temperature for photosynthesis (deg C). See
+##'   Mahadevan et al (2008) eq. 6.
+##' @param Tlow numeric; minimum temperature for respiration (deg C). See
+##'   Mahadevan et al (2008) eq. 10.
+##' @param tower_date chron vector; timestamps for all tower observations
+##'   (NEE_obs, T, PAR)
+##' @param NEE_obs numeric vector; eddy covariance observed net ecosystem
+##'   exchange (NEE, umol m-2 s-1)
 ##' @param T numeric vector; observed air temperature (deg C)
-##' @param PAR numeric vector; observed photosynthetically active
-##' radiation (umol m-2 s-1)
+##' @param PAR numeric vector; observed photosynthetically active radiation
+##'   (umol m-2 s-1)
 ##' @param date_nir chron vector; timestamps for NIR reflectance.
 ##' @param rho_nir numeric vector; NIR reflectance values
-##' @param date_swir  chron vector; timestamps for SWIR reflectance.
-##' @param rho_swir  numeric vector; SWIR reflectance values
-##' @param date_EVI  chron vector; timestamps for enhanced vegetation index (EVI).
+##' @param date_swir chron vector; timestamps for SWIR reflectance.
+##' @param rho_swir numeric vector; SWIR reflectance values
+##' @param date_EVI chron vector; timestamps for enhanced vegetation index
+##'   (EVI).
 ##' @param EVI numeric vector; EVI values.
-##' @param phen factor; phenology dynamics.  levels are ginc (onset
-##' greenness increase), gdec (onset greenness decrease), gmin (onset
-##' greenness minimum), gmax (onset greenness maximum).  If not
-##' specified, phenology dynsamics are calculated from EVI using a
-##' method similar to Zhang et al (2003).
-##' @return an object of class VPRM_driver_data.  Has fields (see
-##' above arguments for definitions): name_long, name_short, lat, lon,
-##' PFT, note, Tmin, Tmax, Topt, Tlow.  The data for NEE_obs, T, PAR,
-##' rho_nir, rho_swir, EVI, and phen should be accessed using the
-##' as.data.frame method.
-##' @references Mahadevan, P., S. C. Wofsy, D. M. Matross, X. Xiao,
-##' A. L. Dunn, J. C. Lin, C. Gerbig, J. W. Munger, V. Y. Chow, and
-##' E. W. Gottlieb (2008), A satellite-based biosphere
-##' parameterization for net ecosystem CO2 exchange: Vegetation
-##' Photosynthesis and Respiration Model (VPRM), Global
-##' Biogeochem. Cycles, 22, GB2005, doi:10.1029/2006GB002735.
-##' @references Xiaoyang Zhang, Mark A. Friedl, Crystal B. Schaaf,
-##' Alan H. Strahler, John C.F. Hodges, Feng Gao, Bradley C. Reed,
-##' Alfredo Huete, Monitoring vegetation phenology using MODIS, Remote
-##' Sensing of Environment, Volume 84, Issue 3, March 2003, Pages
-##' 471-475, ISSN 0034-4257,
-##' http://dx.doi.org/10.1016/S0034-4257(02)00135-9.
+##' @param refEVI numeric vector; reference EVI values. Only required for
+##'   urbanVPRM. See Hardiman et al SI section S2.4.
+##' @param ISA numeric vector; impervious surface area values. Only required for
+##'   urbanVPRM. Must vary between 0.0 and 1.0.
+##' @param LSWI numeric vector; Land Surface Water Index. If not provided will
+##'   be calculated from rho_nir and rho_swir.
+##' @param phen factor; phenology dynamics. levels are ginc (onset greenness
+##'   increase), gdec (onset greenness decrease), gmin (onset greenness
+##'   minimum), gmax (onset greenness maximum). If not specified, phenology
+##'   dynsamics are calculated from EVI using a method similar to Zhang et al
+##'   (2003).
+##' @param model_form string, optional; form of VPRM model to use. Options are
+##'   "Mahadevan07" (default) to use the VPRM formulation of Mahadevan et al.
+##'   (2007), or "urban" to use the urbanVPRM formulation of Hardiman et al.
+##'   (2017). If set to "urban".
+##' @return an object of class VPRM_driver_data. Has fields (see above arguments
+##'   for definitions): name_long, name_short, lat, lon, PFT, note, Tmin, Tmax,
+##'   Topt, Tlow. The data for NEE_obs, T, PAR, rho_nir, rho_swir, EVI, and phen
+##'   should be accessed using the as.data.frame method.
+##' @references Mahadevan, P., S. C. Wofsy, D. M. Matross, X. Xiao, A. L. Dunn,
+##'   J. C. Lin, C. Gerbig, J. W. Munger, V. Y. Chow, and E. W. Gottlieb (2008),
+##'   A satellite-based biosphere parameterization for net ecosystem CO2
+##'   exchange: Vegetation Photosynthesis and Respiration Model (VPRM), Global
+##'   Biogeochem. Cycles, 22, GB2005, doi:10.1029/2006GB002735.
+##' @references Xiaoyang Zhang, Mark A. Friedl, Crystal B. Schaaf, Alan H.
+##'   Strahler, John C.F. Hodges, Feng Gao, Bradley C. Reed, Alfredo Huete,
+##'   Monitoring vegetation phenology using MODIS, Remote Sensing of
+##'   Environment, Volume 84, Issue 3, March 2003, Pages 471-475, ISSN
+##'   0034-4257, http://dx.doi.org/10.1016/S0034-4257(02)00135-9.
 ##' @author Timothy W. Hilton
 ##' @import chron
 ##' @export
@@ -82,20 +91,16 @@
 ##'                            EVI=PFa_evi[['evi']],
 ##'                            phen=NA)
 ##' print(head(as.data.frame(pfa_dd)))
-
-VPRM_driver_data <- function( name_long="",
+VPRM_driver_data <- function(name_long="",
                         name_short="",
                         lat=NA,
                         lon=NA,
                         PFT=NA,
                         note="",
-                        ## scalars
-                        Tmin=0, #minimum respiration temperature (C)
-                        Tmax=40, #maximum respiration temperature (C)
-                        Topt=20, #optimal respiration temperature (C)
+                        Tmin=0,
+                        Tmax=40,
+                        Topt=20,
                         Tlow=2,
-                        ## ---
-                        ## observed fields
                         tower_date=NA,
                         NEE_obs=NA,
                         T=NA,
@@ -106,50 +111,60 @@ VPRM_driver_data <- function( name_long="",
                         rho_swir=NA,
                         date_EVI,
                         EVI=NA,
-                        ## phen should be a data frame with columns "date" and "phen"
-                        phen=NA ) {
+                        refEVI=NA,
+                        ISA=NA,
+                        LSWI=NA,
+                        phen=NA,
+                        model_form='Mahadevan07') {
 
-obj <- list( name_long=name_long,
+## all the 'scalar' fields; that is, fields that are single values, not time
+## series
+obj <- list(name_long=name_long,
             name_short=name_short,
             lat=lat,
             lon=lon,
             note=note,
             PFT=PFT,
             Tmin=Tmin,
             Tmax=Tmax,
-            Topt=Topt )
+            Topt=Topt,
+            Tlow=Tlow,
+            model_form=model_form)
 
 ## interpolate reflectances and EVI onto tower timestamps
-EVI <- interpMODIS( date_EVI, EVI, tower_date, "linear" )[['val']]
-rho_nir <- interpMODIS( date_nir, rho_nir, tower_date, "linear" )[['val']]
-rho_swir <- interpMODIS( date_swir, rho_swir, tower_date, "linear" )[['val']]
+EVI <- interpMODIS(date_EVI, EVI, tower_date, "linear")[['val']]
+rho_nir <- interpMODIS(date_nir, rho_nir, tower_date, "linear")[['val']]
+rho_swir <- interpMODIS(date_swir, rho_swir, tower_date, "linear")[['val']]
 ## if phen contains only NA, then calculate phenology transition dates
 ## from EVI
-if ( identical( all( is.na( phen ) ), TRUE ) ) {
+if (identical(all(is.na(phen)), TRUE)) {
     EVI_df <- data.frame(sitecode=name_short,
                          t=tower_date,
                          EVI=EVI)
-    phen <- detect_large_greenness_change_periods( EVI_df)
+    phen <- detect_large_greenness_change_periods(EVI_df)
 }
 ## interpolate phenology status (greenness increasing, greenness max,
 ## greenness decreasing, greenness minimum) onto tower dates
-phen <- interp_phenology( phen, tower_date )
+phen <- interp_phenology(phen, tower_date)
 
 ## try/catch here!
-obj[['data']] <- data.frame( date=tower_date,
+obj[['data']] <- data.frame(date=tower_date,
                        NEE_obs=NEE_obs,
                        T=T,
                        PAR=PAR,
                        rho_nir=rho_nir,
                        rho_swir=rho_swir,
                        EVI=EVI,
-                       phen=phen[['phen']] )
+                       refEVI=refEVI,
+                       ISA=ISA,
+                       LSWI=LSWI,
+                       phen=phen[['phen']])
 
-class( obj ) <- c( 'VPRM_driver_data', class( obj ) )
+class(obj) <- c('VPRM_driver_data', class(obj))
 
-obj <- calculate_VPRM_derived_input_fields( obj )
+obj <- calculate_VPRM_derived_input_fields(obj)
 
-return( obj )
+return(obj)
 
 }
 
@@ -164,37 +179,43 @@ return( obj )
 ## @return VPRM_driver_data object; same as input argument obj, but with
 ## derived fields filled in.
 ## @author Timothy W. Hilton
-calculate_VPRM_derived_input_fields <- function( obj ) {
+calculate_VPRM_derived_input_fields <- function(obj) {
 
   ## make sure class of obj is correct
-  if ( !( 'VPRM_driver_data' %in% class( obj ) ) )
-    stop( 'obj must be a VPRM_driver_data object' )
-
-  obj[['data']][['LSWI']] <- with( obj[['data']],  getLSWI( rho_nir,
-                                                           rho_swir ) )
-  obj[['data']][['Tscale']] <- getTscale( obj[['data']][['T']],
+  if (!('VPRM_driver_data' %in% class(obj)))
+    stop('obj must be a VPRM_driver_data object')
+
+
+  if (identical(all(is.na(obj[['data']][['LSWI']])), TRUE)) {
+    obj[['data']][['LSWI']] <- with(obj[['data']],
+                                    getLSWI(rho_nir, rho_swir))
+  }
+  obj[['data']][['Tscale']] <- getTscale(obj[['data']][['T']],
                                          obj[['Tmax']],
                                          obj[['Tmin']],
-                                         obj[['Topt']] )
-  if ( any( !is.na( obj[['data']][['phen']] ) ) ) {
-    ## for non-evergreen classes, derive Pscale from phenology
-    obj[['data']][['Pscale']] <- with( obj[['data']], getPscale( LSWI, phen ) )
+                                         obj[['Topt']])
+  if (obj[['model_form']] == 'urban') {
+    obj[['data']][['Pscale']] <- with(obj[['data']], getPscale_urban(EVI))
   } else {
-    ## for evergreen classes, set Pscale to 1.0 (Mahadevan et al, 2008)
-    obj[['data']][['Pscale']] <- 1.0
+    if (any(!is.na(obj[['data']][['phen']]))) {
+      ## for non-evergreen classes, derive Pscale from phenology
+      obj[['data']][['Pscale']] <- with(obj[['data']], getPscale(LSWI, phen))
+    } else {
+      ## for evergreen classes, set Pscale to 1.0 (Mahadevan et al, 2008)
+      obj[['data']][['Pscale']] <- 1.0
+    }
   }
 
-  obj[['data']][['Wscale']] <- with( obj[['data']],
-                                    getWscale( LSWI, max( LSWI, na.rm=TRUE ) ) )
+  obj[['data']][['Wscale']] <- with(obj[['data']],
+                                    getWscale(LSWI, max(LSWI, na.rm=TRUE)))
 
   Tresp <- obj[['data']][['T']]
   Tresp[ Tresp < obj[['Tlow']] ] <- obj[['Tlow']]
   obj[['data']][['Tresp']] <- Tresp
 
-  return( obj )
+  return(obj)
 }
 
-
 ##' method for VPRM_driver_data
 ##'
 ##' The data field of the VPRM_driver_data object is kept.  Other
@@ -226,7 +247,25 @@ calculate_VPRM_derived_input_fields <- function( obj ) {
 ##'                            EVI=PFa_evi[['evi']],
 ##'                            phen=NA)
 ##' print(head(as.data.frame(pfa_dd)))
-as.data.frame.VPRM_driver_data <- function( x, ... ) {
-  return( x[[ 'data' ]] )
+as.data.frame.VPRM_driver_data <- function(x, ...) {
+  return(x[[ 'data' ]])
 }
 
+
+##' calculate photosynthetically available radiation (PAR) from incident
+##' downward shortwave radiation (SWDN)
+##'
+##' uses the mean result from Britton and Dodd (1976) table I. Averages the
+##' tables two values for Apr-May together, then averages all bimonthly
+##' averages.
+##' @title downward shortwave radiation to PAR
+##' @param SWDN numeric vector; downward shortwave radiation (W m-2)
+##' @return numeric vector containing PAR values (umol m-2 s-1)
+##' @author Timothy W. Hilton
+##' @references Britton, C. M. and Dodd, J. D., 1976. Relationships of
+##'   photosynthetically active radiation and shortwave irradiance. Agric.
+##'   Meteorol., 17: 1--7. https://doi.org/10.1016/0002-1571(76)90080-7
+##' @export
+SWDN_2_PAR <- function(SWDN) {
+  return(SWDN * 2.17305)
+}