fix ALB and EMIS approach, clean-up

.gitignore

@@ -0,0 +1 @@
+venv/

README.md

@@ -2,8 +2,11 @@
 
 Set of tools to use Local Climate Zone (LCZ)-based urban canopy parameters in DWD's COSMO-CLM NWP and regional climate model.
 
+## Citaton
+Varentsov, M., Samsonov, T., Demuzere, M., (under review). Impact of urban canopy parameters on a megacity’s modelled thermal environment. Atmosphere.
+
 ## Context
-TERRA_URB, the urban canopy parameterization in COSMO-CLM ([Wouters et al., 2016](https://gmd.copernicus.org/articles/9/3027/2016/)), by default uses impervious surface area information from the [Copernicus Land Monitoring Service](https://land.copernicus.eu/pan-european/high-resolution-layers/imperviousness) (for Europe) / [National Geophysical Data Center](https://databasin.org/datasets/016d2235a5ed43ad83ceeed6c408d149) (global) and anthropogenic heat flux information from [Flanner et al. (2010)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2008gl036465). All other geometrical, thermal and radiative urban canopy parameters are spatially invariant, and set to the values provided in Table 1 of [Wouters et al., 2016](https://gmd.copernicus.org/articles/9/3027/2016/)).
+TERRA_URB is the urban canopy parameterization embedded in TERRA-ML, the land surface model in COSMO-CLM. By default it uses impervious surface area information from the [Copernicus Land Monitoring Service](https://land.copernicus.eu/pan-european/high-resolution-layers/imperviousness) (for Europe) / [National Geophysical Data Center](https://databasin.org/datasets/016d2235a5ed43ad83ceeed6c408d149) (global) and anthropogenic heat flux information from [Flanner et al. (2010)](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2008gl036465). All other geometrical, thermal and radiative urban canopy parameters are spatially invariant, and set to the bulk values provided in Table 1 of [Wouters et al., 2016](https://gmd.copernicus.org/articles/9/3027/2016/)).
 
 The set of tools provided in this repo allow one to introduce LCZ-based urban canopy parameters, compiled from [Stewart and Oke (2012)](http://10.1175/BAMS-D-11-00019.1) and [Stewart et al. (2014)](http://10.1002/joc.3746).
 
@@ -12,30 +15,35 @@ This work is an outcome of AEVUS I and II, the COSMO Priority Tasks on "Analysis
 
 ## Requirements
 * Be a member of the [COSMO-CLM community](https://wiki.coast.hzg.de/clmcom/), in order to be able to access [EXTPAR](https://wiki.coast.hzg.de/clmcom/external-data-98599196.html).
-* have your domain file available from EXTPAR (netcdf file)
-* an LCZ map covering the same region of interest. If not available, please contact me.
-
-
-## Procedure (more info to come)
+* Have your domain file available from EXTPAR (netcdf file)
+* an LCZ map covering the same region of interest. Sources for existing LCZ maps:
+    * Europe: [paper](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0214474) | [data](http://www.wudapt.org/continental-lcz-maps/)
+    * Continental United States: [paper](https://doi.org/10.1038/s41597-020-00605-z) | [data](https://doi.org/10.6084/m9.figshare.11416950.v1) 
+    * An online LCZ Generator tool is currently under development; a beta version can be accessed [here](https://lcz-generator.geographie.rub.de/). Please contact Matthias.Demuzere @ rub.de for more information.
+
 
-`from utils import *`
 
-### Step 1: Assign UCP values to LCZ map and convert to COSMO Grid
-`lcz_to_cosmo(ucpFile, clmFile, lczFile, bandNr, ucpVersion, nrLcz=17,
-              interpMethod='linear', aggregation=True, aggregationScale=2,
-              isaWeight=True, saiWeight=False, fileNameExt='')`
+## Instructions
 
-### Step 2: Address the double counting issue.
-`remove_double_counting(clmFile,gcFile)`
+It is advised to use a python virtual environment:
+1. Go into scriptdir: `cd /SCRIPT/DIR/`
+2. Create virtual environment: `python3 -m venv venv` or `virtualenv venv`
+3. Install module requirements: `venv/bin/pip install -r requirements.txt`
+4. Use `venv/bin/pip/python` to run scripts.
 
+The `requirements.txt` can be generated using `pipreqs`: 
+```
+cd /SCRIPT/DIR/
+pipreqs --ignore=terra/ .
+```
 
 
-## Python version, required packages
-This code was developed and tested using python 3.6.10.
+### Execute
 
-* pandas
-* numpy
-* xarray
-* scipy
+The run code is currently configured for the Moscow case, as developed in Varentsov et al.
+CLM and LCZ input data used in this study is provided under `data/`.
 
+```
+venv/bin/pip/python run.py
+```
 

requirements.txt

@@ -0,0 +1,5 @@
+rasterio==1.1.4
+scipy==1.4.1
+numpy==1.19.4
+xarray==0.15.1
+pandas==1.0.5

run.py

@@ -0,0 +1,36 @@
+import os
+abspath = os.path.abspath(__file__)
+WORKDIR = os.path.dirname(abspath)
+os.chdir(WORKDIR)
+from utils import lcz_to_cosmo, remove_double_counting
+
+## Adjust directory here
+CLMFILE = f"{WORKDIR}/data/MSK_0.009bg3_Globcover.nc"
+LCZFILE = f"{WORKDIR}/data/LCZ_Russia_Moscow.tif"
+
+# FILES
+UCPFILE = f"{WORKDIR}/tables/LCZ_UCP_default.csv"
+GCFILE  = f"{WORKDIR}/tables/globcover_lookup.csv"
+
+
+# EXECUTE FUNCTIONS
+# 1. Assign UCP values to LCZ map and convert to COSMO Grid
+CLM_FILE_NEW = lcz_to_cosmo(
+    ucpFile=UCPFILE,
+    clmFile=CLMFILE,
+    lczFile=LCZFILE,
+    bandNr=3,
+    ucpVersion='default',
+    nrLcz=17,
+    interpMethod='linear',
+    aggregation=True,
+    aggregationScale=2,
+    isaWeight=True,
+    saiWeight=False,
+    fileNameExt='_Varentsov_etal_Atm')
+
+# 2. Address the double counting issue.
+remove_double_counting(
+    clmFile=CLM_FILE_NEW,
+    gcFile=GCFILE,
+    removeUrban=True)

utils.py

@@ -3,6 +3,7 @@
 import xarray as xr
 from scipy.interpolate import RegularGridInterpolator
 import scipy.ndimage
+import rasterio
 
 ## Helper function to prepare the urban canopy data.
 def prepare_ucp_lookup(ucpFile, saiWeight=False, snow_f=0, alb_snow=0.70, emi_snow=0.997):
@@ -11,10 +12,6 @@ def prepare_ucp_lookup(ucpFile, saiWeight=False, snow_f=0, alb_snow=0.70, emi_sn
 
         AUTHOR: Matthias Demuzere (matthias.demuzere [@] rub.de)
 
-        VERSION: 1.0
-
-        DATE: 2020-07-01
-
         :param ucpFile  : Absolute path to urban canopy parameter csv file.
         :param saiWeight: Weigh parameters according to Surface Area Index (Default = False)
         :param snow_f   : snow fraction (default = 0)
@@ -46,57 +43,79 @@ def prepare_ucp_lookup(ucpFile, saiWeight=False, snow_f=0, alb_snow=0.70, emi_sn
     ## Read look-up table
     ucp = pd.read_csv(ucpFile,sep=';',index_col=0).iloc[:17,:]
 
-    ## Canyon albedo reduction factor, eq. 15 Wouters et al. (2016)
+    ## Canyon albedo reduction factor, eq. 15
     psi_canyon = np.exp(-0.6 * ucp['URB_H2W'])
-    psi_canyon[10:] = 0 # Set to zero for non-urban LCZs
+    psi_canyon[10:] = 0  # Set to zero for non-urban LCZs
 
-    ## Total albedo reduction factor, eq. 14 Wouters et al. (2016)
-    #psi_bulk = ucp['URB_BLDFR'] + (1-ucp['URB_BLDFR'])*psi_canyon*ucp['URB_H2W']
+    ## Total albedo reduction factor, eq. 14
+    psi_bulk = psi_canyon * (1 - ucp['URB_BLDFR']) + ucp['URB_BLDFR']
+    psi_bulk[10:] = 0  # Set to zero for non-urban LCZs
 
-    ## Bulk shortwave albedo
+    ## Bulk shortwave albedo, using facet information. Eq 16
     alb_roof_snow = ucp['URB_RfALB'] * (1. - snow_f) + alb_snow * snow_f
     alb_road_snow = ucp['URB_RdALB'] * (1. - snow_f) + alb_snow * snow_f
     alb_wall_snow = ucp['URB_WaALB'] * (1. - snow_f) + alb_snow * snow_f
-    ucp['URB_SALB'] = (alb_road_snow + 2. * ucp['URB_H2W'] * alb_wall_snow) / \
-                            (1. + 2. * ucp['URB_H2W']) * psi_canyon * (1. - ucp['URB_BLDFR']) \
-                            + alb_roof_snow * ucp['URB_BLDFR']
-    ucp.loc[11:,'URB_SALB'] = 0
-    #ucp['URB_TALB'] = ucp['URB_SALB'].copy()
+    ucp['URB_SALB_BK'] = (alb_road_snow + 2. * ucp['URB_H2W'] * alb_wall_snow) / \
+                         (1. + 2. * ucp['URB_H2W']) * psi_canyon * (1. - ucp['URB_BLDFR']) \
+                         + alb_roof_snow * ucp['URB_BLDFR']
+    ucp.loc[11:, 'URB_SALB_BK'] = 0
+    # ucp['URB_TALB'] = ucp['URB_SALB'].copy()
 
-    ## Bulk emissivity
+    ## Bulk emissivity, using facet information. Eq 16
     emi_roof_snow = (1. - ucp['URB_RfEMI']) * (1. - snow_f) + (1. - emi_snow) * snow_f
     emi_road_snow = (1. - ucp['URB_RdEMI']) * (1. - snow_f) + (1. - emi_snow) * snow_f
     emi_wall_snow = (1. - ucp['URB_WaEMI']) * (1. - snow_f) + (1. - emi_snow) * snow_f
-    ucp['URB_EMIS'] = 1. - ((emi_road_snow + 2. * ucp['URB_H2W'] * emi_wall_snow) \
-                     / (1. + 2. * ucp['URB_H2W']) * psi_canyon * (1. - ucp['URB_BLDFR']) \
-                     + emi_roof_snow * ucp['URB_BLDFR'])
-    ucp.loc[11:,'URB_EMIS'] = 0
+    ucp['URB_EMIS_BK'] = 1. - ((emi_road_snow + 2. * ucp['URB_H2W'] * emi_wall_snow) \
+                               / (1. + 2. * ucp['URB_H2W']) * psi_canyon * (1. - ucp['URB_BLDFR']) \
+                               + emi_roof_snow * ucp['URB_BLDFR'])
+    ucp.loc[11:, 'URB_EMIS_BK'] = 0
 
     ## Bulk thermal albedo
-    ucp['URB_TALB'] = 1- ucp['URB_EMIS']
-    ucp.loc[11:, 'URB_TALB'] = 0
-
+    ucp['URB_TALB_BK'] = 1 - ucp['URB_EMIS_BK']
+    ucp.loc[11:, 'URB_TALB_BK'] = 0
 
     ## Calculate Surface Area Index from geometrical considerations (Eq. 3)
     SAI = (1. + 2. * ucp['URB_H2W']) * (1. - ucp['URB_BLDFR']) + ucp['URB_BLDFR']
 
-    ## Get bulk Heat capacity and conductivity, using eq. 10, 11 and 4.
-    ucp['URB_HCON'] = ((1-ucp['URB_BLDFR']) / SAI ) * \
-              (2*ucp['URB_H2W']*ucp['URB_WaHCON'] + ucp['URB_RdHCON']) + \
-              ( ucp['URB_BLDFR'] / SAI * ucp['URB_RfHCON'])
-    ucp['URB_HCAP'] = ((1-ucp['URB_BLDFR']) / SAI ) * \
-              (2*ucp['URB_H2W']*ucp['URB_WaHCAP'] + ucp['URB_RdHCAP']) + \
-              ( ucp['URB_BLDFR'] / SAI * ucp['URB_RfHCAP'])
+    ## Get mean Heat capacity and conductivity, using eq. 10, 11 and 4.
+    ucp['URB_HCON'] = ((1 - ucp['URB_BLDFR']) / SAI) * \
+                      (2 * ucp['URB_H2W'] * ucp['URB_WaHCON'] + ucp['URB_RdHCON']) + \
+                      (ucp['URB_BLDFR'] / SAI * ucp['URB_RfHCON'])
+    ucp['URB_HCAP'] = ((1 - ucp['URB_BLDFR']) / SAI) * \
+                      (2 * ucp['URB_H2W'] * ucp['URB_WaHCAP'] + ucp['URB_RdHCAP']) + \
+                      (ucp['URB_BLDFR'] / SAI * ucp['URB_RfHCAP'])
+
+    ## Mean facet-level albedo and emissivity based on eq. 10
+    ## Only added for testing and potential comparison with other models
+    ## These values are currently not used in TERRA_URB.
+    ucp['URB_EMIS_FL'] = ((1 - ucp['URB_BLDFR']) / SAI) * \
+                         (2 * ucp['URB_H2W'] * ucp['URB_WaEMI'] + ucp['URB_RdEMI']) + \
+                         (ucp['URB_BLDFR'] / SAI * ucp['URB_RfEMI'])
+    ucp['URB_SALB_FL'] = ((1 - ucp['URB_BLDFR']) / SAI) * \
+                         (2 * ucp['URB_H2W'] * ucp['URB_WaALB'] + ucp['URB_RdALB']) + \
+                         (ucp['URB_BLDFR'] / SAI * ucp['URB_RfALB'])
+    ucp['URB_TALB_FL'] = 1 - ucp['URB_EMIS_FL']
+    ucp.loc[11:, 'URB_SALB_FL'] = 0
+    ucp.loc[11:, 'URB_TALB_FL'] = 0
+
+    ## For now, TERRA-URB only reads in one average facet-level albedo.
+    ## The bulk calculation from eqs. 13 is done within TERRA_URB
+    ## Therefore, the bulk value needs to be reversed back to a mean
+    ## facet value, so that eq. 13 is solved for alb = alb_bulk / psi_bulk
+    ## The same is done for the emissivity.
+    ucp['URB_SALB'] = ucp['URB_SALB_BK'] / psi_bulk
+    ucp['URB_TALB'] = ucp['URB_TALB_BK'] / psi_bulk
+    ucp['URB_EMIS'] = 1 - ucp['URB_TALB']
 
     ## Also add the thermal admittance
-    #ucp['URB_TADM'] = (ucp['URB_HCAP']*ucp['URB_HCON'])**0.5
+    # ucp['URB_TADM'] = (ucp['URB_HCAP']*ucp['URB_HCON'])**0.5
 
     ## iS SAI weighting requested, according to Eq. 4?
     ## This is done within TERRA_URB, so no need to do for COSMO/CLM input files.
     if saiWeight:
         ucp['URB_HCON'] = ucp['URB_HCON'] * SAI
         ucp['URB_HCAP'] = ucp['URB_HCAP'] * SAI
-        #ucp['URB_TADM'] = ucp['URB_TADM'] * SAI
+        # ucp['URB_TADM'] = ucp['URB_TADM'] * SAI
 
     return ucp
 
@@ -108,10 +127,6 @@ def cosmo_interpolator(xLcz, yLcz, dataLcz, xClm, yClm, interpMethod='linear',
     """
         AUTHOR: Matthias Demuzere (matthias.demuzere [@] rub.de)
 
-        VERSION: 1.0
-
-        DATE: 2020-07-01
-
         :param xLcz: values of LCZ longitudes
         :param yLcz: values of LCZ latitudes
         :param dataLcz: 2D array of LCZ parameter value
@@ -151,10 +166,6 @@ def lcz_to_cosmo(ucpFile, clmFile, lczFile, bandNr, ucpVersion, nrLcz=17,
 
         AUTHOR: Matthias Demuzere (matthias.demuzere [@] rub.de)
 
-        VERSION: 1.0
-
-        DATE: 2020-07-01
-
         :param ucpFile: full absolute path name to ucp .csv table.
         :param clmFile: full absolute path name to COSMO-CLM domain file
         :param lczFile: full absolute path name to lcz geotiff file.
@@ -217,6 +228,12 @@ def lcz_to_cosmo(ucpFile, clmFile, lczFile, bandNr, ucpVersion, nrLcz=17,
                      'URB_SALB',
                      'URB_TALB',
                      'URB_EMIS',
+                     'URB_SALB_FL',
+                     'URB_TALB_FL',
+                     'URB_EMIS_FL',
+                     'URB_SALB_BK',
+                     'URB_TALB_BK',
+                     'URB_EMIS_BK',
                      'URB_HCON',
                      'URB_HCAP']
 
@@ -294,10 +311,6 @@ def remove_double_counting(clmFile,gcFile,removeUrban=True,qLow=0.25,qHigh=0.75,
 
         AUTHOR: Matthias Demuzere (matthias.demuzere [@] rub.de)
 
-        VERSION: 1.0
-
-        DATE: 2020-07-01
-
         :param clmFile: full absolute path name to COSMO-CLM domain file
         :param gcFile: full absolute path name to globcover look-up table
         :param removeUrban: Boolean to indicate if URBAN effect from GlobCover needs to be removed (default=True)