From c2972504020438ff8a26bf4f2240bf65b81f0b44 Mon Sep 17 00:00:00 2001 From: ted Date: Tue, 23 Apr 2019 22:17:57 -0700 Subject: [PATCH] Added docs and tools for procedure to generate projected land cover change parameters. --- README.md | 4 +- docs/Procedure.md | 71 +-------- docs/Procedure_Projected_Land_Cover_Change.md | 46 ++++++ docs/Utilities.md | 69 ++++++++ tools/aggregate_lc_map.py | 147 ++++++++++++++++++ tools/gapfill_vic_params.py | 34 ++-- 6 files changed, 283 insertions(+), 88 deletions(-) create mode 100644 docs/Procedure_Projected_Land_Cover_Change.md create mode 100644 docs/Utilities.md create mode 100755 tools/aggregate_lc_map.py diff --git a/README.md b/README.md index 1605024..e4d31f1 100644 --- a/README.md +++ b/README.md @@ -51,7 +51,9 @@ These parameter files were designed for use with VIC 5 (image driver). VIC 5 ima - README.md - this file docs/ - - [Procedure.md](./docs/Procedure.md) - list of processing steps, the scripts that run them, and their usage + - [Procedure.md](./docs/Procedure.md) - Main Procedure; list of processing steps, the scripts that run them, and their usage + - [Procedure_Projected_Land_Cover_Change.md](./docs/Procedure_Projected_Land_Cover_Change.md) - Procedure for modifying existing VIC parameters to account for future projected land cover change; list of processing steps, the scripts that run them, and their usage + - [Utilities.md](./docs/Utilities.md) - Descriptions of utility scripts, and their usage examples/ - batch files containing example commands and arguments diff --git a/docs/Procedure.md b/docs/Procedure.md index 68fe172..8a7d0ec 100644 --- a/docs/Procedure.md +++ b/docs/Procedure.md @@ -158,7 +158,7 @@ Before stage 4 can begin, we need 2 other files: For the CONUS_MX and USMX domains, domain files are available for download on Zenodo [(PITRI Precipitation Disaggregation Parameters)](https://zenodo.org/record/2564019). For these same domains, "source" parameter files based on Livneh et al (2015) ("L2015" henceforth) are availabled for download on Zenodo [(MOD-LSP VIC Parameters)](https://zenodo.org/record/2612560). -If you have a different domain, you can either clip out your domain from within these domain and parameter files (if your domain is inside those domains) or create a domain file from elemental inputs (mask, and DEM) with a utility script (see the "Utility Scripts" section below). To create a NetCDF VIC-5 compliant "source" parameter file for a domain outside L2015, you will need ascii-format VIC soil, snow, and vegetation parameters (and vegetation library file) for the domain, and you can use the tool [tonic](https://github.com/UW-Hydro/tonic). +If you have a different domain, you can either clip out your domain from within these domain and parameter files (if your domain is inside those domains) or create a domain file from elemental inputs (mask, and DEM) with a utility script (see the [Utilities](Utilities.md) documentation). To create a NetCDF VIC-5 compliant "source" parameter file for a domain outside L2015, you will need ascii-format VIC soil, snow, and vegetation parameters (and vegetation library file) for the domain, and you can use the tool [tonic](https://github.com/UW-Hydro/tonic). ### Option 1: if you have divided your domain into 10x10 degree boxes @@ -284,72 +284,3 @@ This stage is only necessary if you wish to drive VIC with an explicit 17-year t `$PREFIX` = prefix of input/output filenames. `$OUTDIR` = output directory. -## Utility Scripts - -Some useful utility scripts: - - To create a domain file for a new domain: - - `create_domain_file_from_asc_inputs.py -m $MASK -f $FRACMASK -d $DEM -p $PRCP_PARAMDIR -o $OUTFILE` - - where - - `$MASK` = ESRI ascii-format mask delineating the domain (1 = in domain; 0 = outside) - `$FRACMASK` = ESRI ascii-format grid showing the fraction of each cell's area that lies within the domain - `$DEM` = ESRI ascii-format grid of elevation values (Digital Elevation Model) - `$PRCP_PARAMDIR` = directory containing ESRI ascii-format grid files of precipitation parameters (2 parameters, "duration" and "peak_time", with 12 monthly mean values each, one file per month) - `$OUTFILE` = output path/filename - - - To create initial state files for MetSim: - - `create_metsim_state_file_from_domain_file.py -i $INFILE -s $STATEDATE -o $OUTFILE` - - where - - `$INFILE` = domain file - `$STATEDATE` = date of 1 day prior to the desired start; format `$YYYY-$MM-$SDD` - `$OUTFILE` = output MetSim state file - - - To clip out a spatial subset (lat/lon box) of a netcdf file: - - `gridclip.py -i $INFILE -s $SOUTH -n $NORTH -w $WEST -e $EAST -o $OUTFILE` - - where - - `$INFILE` = input file - `$SOUTH`,`$NORTH`,`$WEST`,`$EAST` = geographic boundaries of subset - `$OUTFILE` = output file - - - To subsample VIC parameters to a finer resolution AND/OR clip to a new irregularly-shaped domain (via the mask of a domain file): - - `subsample_vic_params_over_domain.py -i $INFILE -d $DOMAINFILE [-n] [-s] [-f] -o $OUTFILE` - - where - - `$INFILE` = input file - `$DOMAINFILE` = domain file to clip to (at either the same or finer resolution) - `-n` = optional; if specified, use default fill ("nodata") values - `-s` = optional; if specified, `$INFILE` contains snow bands - `-f` = optional; if specified, `$INFILE` contains Fcanopy - `$OUTFILE` = output file - - - To subsample VIC gridded daily forcings to a finer resolution AND/OR clip to a new irregularly-shaped domain (via the mask of a domain file): - - `subsample_forcings_over_domain.py -i $INFILE -d $DOMAINFILE -v $VARNAMELIST -o $OUTFILE` - - where - - `$INFILE` = input file - `$DOMAINFILE` = domain file to clip to (at either the same or finer resolution) - `$VARNAMELIST` = (optional) comma-separated list of variables to process; default = "Prec,Tmin,Tmax,wind" - `$OUTFILE` = output file - - - To set the "run_cell" variable equal to the mask of a forcing file (to ensure that VIC doesn't try to run where no forcings are available): - - `set_run_cell.py -p $PARAMFILE -r $RUNCELLFILE -v $VARNAME -o $OUTFILE` - - where - - `$PARAMFILE` = VIC parameter file - `$RUNCELLFILE` = file from which to get the values for setting "run_cell", e.g., a domain file, or maybe another VIC parameter file - `$VARNAME` = name of variable in `$RUNCELLFILE` the values of which will be assigned to "run_cell" - `$OUTFILE` = output file diff --git a/docs/Procedure_Projected_Land_Cover_Change.md b/docs/Procedure_Projected_Land_Cover_Change.md new file mode 100644 index 0000000..ee5a5fc --- /dev/null +++ b/docs/Procedure_Projected_Land_Cover_Change.md @@ -0,0 +1,46 @@ +# Procedure for Future Projected Land Cover Change + +This document describes the procedure for generating VIC parameters for projected future land cover change. This procedure assumes we have started with a set of land cover maps corresponding to at least one historical map and at least one future projection using the same classification scheme. The historical map should be used to generate historical VIC parameters via the "NLCD_INEGI" version of the [main procedure](Procedure.md). + +The future parameters are then generated by (1) aggregating the future land cover map (ascii-format ESRI grid file) to area fractions at the desired resolution; (2) replacing the `Cv` variable in the historical VIC parameter file with the area fractions from step (1); and (3) filling gaps in the parameters (cells for which classes have appeared in the future that did not exist there in the historical map) by spatially interpolating the parameters for those classes from neighboring cells. + +## 1. Aggregate land cover map + +This step aggregates a land cover map (ascii-format ESRI grid file) to area fractions at the desired resolution, in NetCDF format. The area fractions are stored in a variable named `Cv`. Run the following script: + +`aggregate_lc_map.py -d $DOMAIN_FILE -m $LC_MAP_FILE -t $LC_TABLE -o $OUTFILE` + +where + +`$DOMAIN_FILE` = VIC-5-image-mode-compliant domain file (NetCDF) +`$LC_MAP_FILE` = ESRI ascii-format land cover map +`$LC_TABLE` = Ascii comma-separated-value (csv) file listing the land cover classes, and their names, similar to [../data/USMX/lc_table.USMX.NLCD_INEGI.csv](../data/USMX/lc_table.USMX.NLCD_INEGI.csv) +`$OUTFILE` = output path/filename + +## 2. Replace Cv in historical parameter file + +This step replaces the `Cv` variable of a historical VIC parameter file with the `Cv` variable generated in step 1. Run the following script: + +`replace_cv.py -i $INFILE -c $CV_FILE -o $OUTFILE` + +where + +`$INFILE` = VIC-5-image-mode-compliant historical parameter file +`$CV_FILE` = NetCDF file containing new `Cv` variable, created in step 1 +`$OUTFILE` = output path/filename (this will be the historical VIC parameters but with the future `Cv`) + +## 3. Fill gaps + +This step fills gaps in the parameters (arising from classes appearing in cells where no parameters exist for them) by spatial interpolation from neghboring cells. Run the following script: + +`gapfill_vic_params.py -i $INFILE [-a] [-r] [-m] -o $OUTFILE` + +where + +`$INFILE` = modified VIC parameter file created in step 2 +`-a` = optional flag specifying that a `max_snow_albedo` variable is present +`-r` = optional flag specifying that irrigation variables are present +`-m` = optional flag specifying that impervious area fraction variables are present +`$OUTFILE` = output path/filename for gap-filled parameters + +NOTE: the `-r` and `-m` options are only applicable to the `feature/irrig.imperv.deep_esoil` branch of the `tbohn/VIC` fork of the `UW-Hydro/VIC` GitHub repo, available at [https://github.com/tbohn/VIC/tree/feature/irrig.imperv.deep_esoil](https://github.com/tbohn/VIC/tree/feature/irrig.imperv.deep_esoil). diff --git a/docs/Utilities.md b/docs/Utilities.md new file mode 100644 index 0000000..36cd632 --- /dev/null +++ b/docs/Utilities.md @@ -0,0 +1,69 @@ +# Utility Scripts for the VIC_Landcover_MODIS_NLCD_INEGI Project + +Some useful utility scripts: + - To create a domain file for a new domain: + + `create_domain_file_from_asc_inputs.py -m $MASK -f $FRACMASK -d $DEM -p $PRCP_PARAMDIR -o $OUTFILE` + + where + + `$MASK` = ESRI ascii-format mask delineating the domain (1 = in domain; 0 = outside) + `$FRACMASK` = ESRI ascii-format grid showing the fraction of each cell's area that lies within the domain + `$DEM` = ESRI ascii-format grid of elevation values (Digital Elevation Model) + `$PRCP_PARAMDIR` = directory containing ESRI ascii-format grid files of precipitation parameters (2 parameters, "duration" and "peak_time", with 12 monthly mean values each, one file per month) + `$OUTFILE` = output path/filename + + - To create initial state files for MetSim: + + `create_metsim_state_file_from_domain_file.py -i $INFILE -s $STATEDATE -o $OUTFILE` + + where + + `$INFILE` = domain file + `$STATEDATE` = date of 1 day prior to the desired start; format `$YYYY-$MM-$SDD` + `$OUTFILE` = output MetSim state file + + - To clip out a spatial subset (lat/lon box) of a netcdf file: + + `gridclip.py -i $INFILE -s $SOUTH -n $NORTH -w $WEST -e $EAST -o $OUTFILE` + + where + + `$INFILE` = input file + `$SOUTH`,`$NORTH`,`$WEST`,`$EAST` = geographic boundaries of subset + `$OUTFILE` = output file + + - To subsample VIC parameters to a finer resolution AND/OR clip to a new irregularly-shaped domain (via the mask of a domain file): + + `subsample_vic_params_over_domain.py -i $INFILE -d $DOMAINFILE [-n] [-s] [-f] -o $OUTFILE` + + where + + `$INFILE` = input file + `$DOMAINFILE` = domain file to clip to (at either the same or finer resolution) + `-n` = optional; if specified, use default fill ("nodata") values + `-s` = optional; if specified, `$INFILE` contains snow bands + `-f` = optional; if specified, `$INFILE` contains Fcanopy + `$OUTFILE` = output file + + - To subsample VIC gridded daily forcings to a finer resolution AND/OR clip to a new irregularly-shaped domain (via the mask of a domain file): + + `subsample_forcings_over_domain.py -i $INFILE -d $DOMAINFILE -v $VARNAMELIST -o $OUTFILE` + + where + + `$INFILE` = input file + `$DOMAINFILE` = domain file to clip to (at either the same or finer resolution) + `$VARNAMELIST` = (optional) comma-separated list of variables to process; default = "Prec,Tmin,Tmax,wind" + `$OUTFILE` = output file + + - To set the "run_cell" variable equal to the mask of a forcing file (to ensure that VIC doesn't try to run where no forcings are available): + + `set_run_cell.py -p $PARAMFILE -r $RUNCELLFILE -v $VARNAME -o $OUTFILE` + + where + + `$PARAMFILE` = VIC parameter file + `$RUNCELLFILE` = file from which to get the values for setting "run_cell", e.g., a domain file, or maybe another VIC parameter file + `$VARNAME` = name of variable in `$RUNCELLFILE` the values of which will be assigned to "run_cell" + `$OUTFILE` = output file diff --git a/tools/aggregate_lc_map.py b/tools/aggregate_lc_map.py new file mode 100755 index 0000000..1f17fea --- /dev/null +++ b/tools/aggregate_lc_map.py @@ -0,0 +1,147 @@ +#!/usr/bin/env python + +import os.path +import numpy as np +import xarray as xr +import sys, getopt + +def read_ascii_gridfile(filename, dtype): + + # Read header + f = open(filename, 'r') + header1 = f.readline() + header2 = f.readline() + header3 = f.readline() + header4 = f.readline() + header5 = f.readline() + header6 = f.readline() + ncols = int(header1.partition(' ')[2]) + nrows = int(header2.partition(' ')[2]) + xllcorner = float(header3.partition(' ')[2]) + yllcorner = float(header4.partition(' ')[2]) + cellsize = float(header5.partition(' ')[2]) + nodata = float(header6.partition(' ')[2]) + f.close() + + # Allocate data + data = np.empty((nrows,ncols), dtype=dtype) + + # Read data + tmp = np.loadtxt(filename,dtype=dtype,delimiter =' ', skiprows=6) + for i in range(nrows): +# data[nrows-1-i] = tmp[i].copy() + data[i] = tmp[i].copy() + + llcorner = [yllcorner, xllcorner] + + return (data, nrows, ncols, llcorner, cellsize, nodata) + + +def main(): + domainfile = '' + lc_map_file = '' + lc_table = '' + outfile = '' + try: + opts, args = getopt.getopt(sys.argv[1:],"hd:m:t:o:",["domainfile=","lc_map_file=","lc_table=","outfile="]) + except getopt.GetoptError: + print(sys.argv[0], ' -d -m -t -o ') + sys.exit(2) + for opt, arg in opts: + if opt == '-h': + print(sys.argv[0], ' -d -m -t -o ') + sys.exit() + elif opt in ("-d", "--domainfile"): + domainfile = arg + elif opt in ("-m", "--lc_map_file"): + lc_map_file = arg + elif opt in ("-t", "--lc_table"): + lc_table = arg + elif opt in ("-o", "--outfile"): + outfile = arg + + # Read domain file + ds = xr.open_dataset(domainfile) + + nLat = ds['mask'].shape[0] + nLon = ds['mask'].shape[1] + + lat = np.empty([nLat]) + lon = np.empty([nLon]) + mask = np.empty([nLat,nLon]) + + lat[:] = ds['lat'][:] + lon[:] = ds['lon'][:] + mask[:] = ds['mask'][:] + + cellsize = round((np.max(lat) - np.min(lat)) / float(nLat - 1),6) + minlat = round(np.min(lat) - (0.5 * cellsize),6) + minlon = round(np.min(lon) - (0.5 * cellsize),6) + maxlat = round(np.max(lat) + (0.5 * cellsize),6) + maxlon = round(np.max(lon) + (0.5 * cellsize),6) + + # Read LC table + lc_data = np.loadtxt(lc_table, dtype=bytes, delimiter=',').astype(str) + nClass = lc_data.shape[0] + v_of_classID = {} + classID = np.empty([nClass]).astype(int) + for v in range(nClass): + classID[v] = int(lc_data[v,0]) + v_of_classID[classID[v]] = v + + # LC map + (lc_map, nrows_lc, ncols_lc, llcorner_lc, cellsize_lc, nodata_lc) = \ + read_ascii_gridfile(lc_map_file, np.int32) + + # Aggregate lc classes over domain + Cv = np.zeros([nClass,nLat,nLon]) + count_tot = np.zeros([nLat,nLon]) + for i in range(nrows_lc): + lat_lc = llcorner_lc[0] + (i + 0.5) * cellsize_lc + row = int( (lat_lc - minlat) / cellsize) + print(i,lat_lc,row,lat[row]) + if row < 0 or row >= nLat: + continue + for j in range(ncols_lc): + lon_lc = llcorner_lc[1] + (j + 0.5) * cellsize_lc + col = int( (lon_lc - minlon) / cellsize) + if col < 0 or col >= nLon: + continue + if (mask[row,col] == 1 and lc_map[i,j] != nodata_lc and + lc_map[i,j] in v_of_classID.keys()): + v = v_of_classID[lc_map[i,j]] + Cv[v,row,col] += 1 + count_tot[row,col] += 1 + + for v in range(nClass): + Cv[v] /= count_tot + + # Write output file + ds_out = xr.Dataset( + { + 'mask': (['lat','lon'], mask), + 'Cv': (['class','lat','lon'], Cv), + }, + coords = { + 'lat': (['lat'], lat), + 'lon': (['lon'], lon), + 'class': (['class'], classID), + } + ) + + for varname in ['lat','lon','mask']: + ds_out[varname].attrs = ds[varname].attrs + + ds_out['Cv'].attrs['long_name'] = 'Area fraction for each land cover class' + + for varname in ['lat','lon','mask','Cv']: + ds_out[varname].encoding['zlib'] = True + + print('writing to',outfile) + ds_out.to_netcdf(outfile) + + ds_out.close() + ds.close() + +if __name__ == "__main__": + main() diff --git a/tools/gapfill_vic_params.py b/tools/gapfill_vic_params.py index 6942c03..19b586f 100755 --- a/tools/gapfill_vic_params.py +++ b/tools/gapfill_vic_params.py @@ -10,26 +10,26 @@ def main(): infile = '' outfile = '' - has_nsa = False - has_irr = False - has_imp = False + has_max_snow_albedo = False + has_irrigation = False + has_impervious = False try: - opts, args = getopt.getopt(sys.argv[1:],"hi:nrmo:",["infile=","outfile="]) + opts, args = getopt.getopt(sys.argv[1:],"hi:armo:",["infile=","outfile="]) except getopt.GetoptError: - print(sys.argv[0], ' -i [-n] [-r] [-m] -o ') + print(sys.argv[0], ' -i [-a] [-r] [-m] -o ') sys.exit(2) for opt, arg in opts: if opt == '-h': - print(sys.argv[0], ' -i [-n] [-r] [-m] -o ') + print(sys.argv[0], ' -i [-a] [-r] [-m] -o ') sys.exit() elif opt in ("-i", "--infile"): infile = arg - elif opt in ("-n", "--has_nsa"): - has_nsa = True - elif opt in ("-r", "--has_irr"): - has_irr = True - elif opt in ("-m", "--has_imp"): - has_imp = True + elif opt in ("-a", "--has_max_snow_albedo"): + has_max_snow_albedo = True + elif opt in ("-r", "--has_irrigation"): + has_irrigation = True + elif opt in ("-m", "--has_impervious"): + has_impervious = True elif opt in ("-o", "--outfile"): outfile = arg @@ -59,8 +59,8 @@ def main(): 'annual_prec', 'fs_active', ] - if (has_nsa): - varnames_2d_soil.append('new_snow_albedo') + if (has_max_snow_albedo): + varnames_2d_soil.append('max_snow_albedo') varnames_2d_veg = [ 'Nveg', @@ -99,11 +99,11 @@ def main(): 'trunk_ratio', ] - if (has_irr): + if (has_irrigation): varnames_3d_veg.append('irr_active') varnames_3d_veg.append('ithresh') varnames_3d_veg.append('itarget') - if (has_imp): + if (has_impervious): varnames_3d_veg.append('fimp') varnames_3d_veg.append('feffimp') @@ -119,7 +119,7 @@ def main(): 'veg_rough', 'displacement', ] - if (has_irr): + if (has_irrigation): varnames_4d_clim.append('fcrop') varnames_4d_clim.append('firr') varnames_4d_clim.append('irr_clim')