Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Added docs and tools for procedure to generate projected land cover c…
…hange parameters.
- Loading branch information
ted
committed
Apr 24, 2019
1 parent
27a4684
commit c297250
Showing
6 changed files
with
283 additions
and
88 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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). |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 <domainfile> -m <lc_map_file> -t <lc_table> -o <outfile>') | ||
| sys.exit(2) | ||
| for opt, arg in opts: | ||
| if opt == '-h': | ||
| print(sys.argv[0], ' -d <domainfile> -m <lc_map_file> -t <lc_table> -o <outfile>') | ||
| 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() |
Oops, something went wrong.