diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f7275bb --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +venv/ diff --git a/README.md b/README.md index a86700b..4917987 100644 --- a/README.md +++ b/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 +``` diff --git a/__pycache__/utils.cpython-36.pyc b/__pycache__/utils.cpython-36.pyc index f5aa0a0..8c3ef66 100644 Binary files a/__pycache__/utils.cpython-36.pyc and b/__pycache__/utils.cpython-36.pyc differ diff --git a/data/LCZ_Russia_Moscow.tif b/data/LCZ_Russia_Moscow.tif new file mode 100644 index 0000000..52a449a Binary files /dev/null and b/data/LCZ_Russia_Moscow.tif differ diff --git a/data/MSK_0.009bg3_Globcover.nc b/data/MSK_0.009bg3_Globcover.nc new file mode 100644 index 0000000..a8a62f7 Binary files /dev/null and b/data/MSK_0.009bg3_Globcover.nc differ diff --git a/data/MSK_0.009bg3_Globcover_lcz_3_default_Varentsov_etal_Atm.nc b/data/MSK_0.009bg3_Globcover_lcz_3_default_Varentsov_etal_Atm.nc new file mode 100644 index 0000000..77ad4e0 Binary files /dev/null and b/data/MSK_0.009bg3_Globcover_lcz_3_default_Varentsov_etal_Atm.nc differ diff --git a/data/MSK_0.009bg3_Globcover_lcz_3_default_Varentsov_etal_Atm_fixDC_True.nc b/data/MSK_0.009bg3_Globcover_lcz_3_default_Varentsov_etal_Atm_fixDC_True.nc new file mode 100644 index 0000000..ec16c75 Binary files /dev/null and b/data/MSK_0.009bg3_Globcover_lcz_3_default_Varentsov_etal_Atm_fixDC_True.nc differ diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..92496d8 --- /dev/null +++ b/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 diff --git a/run.py b/run.py new file mode 100644 index 0000000..6a1bfc7 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/utils.py b/utils.py index 289d6f2..6e3f99b 100644 --- a/utils.py +++ b/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)