This notebook will rely on the WRF-Hydro GIS Pre-processing tools, found here:
Throughout this exercise, we will use python variables to store directory paths and other variables. However, we will call all GIS Pre-processing functionality as though it were on the command line. This is done by adding !
syntax before the command-line syntax, to execute the line using bash.
In this cell, Python variables are created that point to the file paths of the test-case data and an output directory is defined to store the data created by these tools.
# Import python core modules
import os
import shutil
# Set root directory for GIS lesson
gis_data_folder = "/home/docker/GIS_Training"
# Change the directory to the GIS_Training directory and get current working directory
os.chdir(gis_data_folder)
cwd = os.getcwd()
# Set paths to known input and output directories and files
data_folder = os.path.join(cwd, 'Croton_Lambert')
in_geogrid = os.path.join(data_folder, 'geo_em.d01.nc')
output_folder = os.path.join(cwd, 'Outputs')
# Clear any outputs from previous runs by deleting (if necessary) and re-creating the output directory
if os.path.exists(output_folder):
shutil.rmtree(output_folder)
os.mkdir(output_folder)
The tool Create_Domain_Boundary_Shapefile.py
takes a WRF (WPS) output file, aka "Geogrid file", and creates a polygon shapefile that defines the boundary of the domain as a single rectangular polygon in projected coordinates. The script will read metadata in the geogrid file and the output shapefile will be in the projection of the WRF domain. The unstaggered grid, or "Mass" grid (e.g. "HGT_M" variable), is used as the routing grid domain by WRF-Hydro.
This is an example of the syntax for calling the Create_Domain_Boundary_Shapefile.py
tool on the command line. By following the tool with -h
or --help
, we are able to call the help arguement which explains the purpose of the tool, shows the different arguements we can use for this tool, as well as the descriptions for each arguement.
When the tool is run in the terminal, it is not necessary to use the exclamation point. In Jupyter Notebook, we can execute command-line syntax using "!".
# Execute script on the command-line, requesting tool help (parameter -h)
! python Create_Domain_Boundary_Shapefile.py -h
Downloading WhiteboxTools pre-compiled binary for first time use ... Decompressing WhiteboxTools_linux_amd64.tar.xz ... WhiteboxTools package directory: /home/docker/miniconda3/lib/python3.8/site-packages/whitebox Downloading testdata ... Script initiated at Thu Nov 5 17:59:38 2020 usage: Create_Domain_Boundary_Shapefile.py [-h] -i IN_NC -o OUT_DIR This tool takes an WRF Geogrid file and creates a single polygon shapefile that makes up the boundary of the domain of the M-grid (HGT_M, for example). optional arguments: -h, --help show this help message and exit -i IN_NC Path to WPS geogrid (geo_em.d0*.nc) file or WRF-Hydro Fulldom_hires.nc file. -o OUT_DIR Output directory.
You will see in the messages above that the tool provides a brief explanation of the expected input and output parameters. This tool requires a geogrid file as input (-i
) and a directory to write the outputs into (-o
)
Now that we know what arguements are needed for this tool, we can enter those arguements and run the tool. For this script, we only need to specify the file path to the WPS geogrid file and an output folder to save the result. The result of this tool is a shapefile that shows the geographic boundary of the domain, as defined in the geogrid file.
When running this tool in Jupyter, we can use brackets around our python variable names, and Jupyter will substitute the variable values when executing the syntax. This is akin to using an environment variable on the command-line. For the sake of repeatability, we also print the full syntax for reference. This can be copied into the terminal if desired.
# Print information to screen for reference
print('Command to run:\n')
print('python Create_Domain_Boundary_Shapefile.py \\\n\t -i {0} \\\n\t -o {1}\n'.format(in_geogrid, output_folder))
# Run the script with required parameters
! python Create_Domain_Boundary_Shapefile.py -i {in_geogrid} -o {output_folder}
Command to run: python Create_Domain_Boundary_Shapefile.py \ -i /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/geo_em.d01.nc \ -o /home/docker/.mnt/common/mollymca/GIS_Training/Outputs Script initiated at Thu Nov 5 18:02:51 2020 WPS netCDF projection identification initiated... Map Projection: Lambert Conformal Conic Using MOAD_CEN_LAT for latitude of origin. Using Standard Parallel 2 in Lambert Conformal Conic map projection. Geo-referencing step completed without error in 0.05 seconds. Created projection definition from input NetCDF GEOGRID file. ESRI Shapefile driver is available. Done producing output vector polygon shapefile in 0.12 seconds Output shapefile: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/geo_em.d01_boundary.shp Process completed in 0.34 seconds.
The messages returned by the tool can be quite useful. You will see the coordinate system information printed to the screen and any other progress messages.
Now that the domain boundary shapefile has been created, we want to see where the domain is relative to other features on a map. The next cell creates a map and adds the domain boundary as a layer. Use the map to explore the domain. A swipe feature allows the basemap to be changed between OpenStreetMap and satellite imagery.
import json
import geopandas
from ipyleaflet import Map, GeoJSON, ScaleControl, FullScreenControl, basemaps, SplitMapControl, basemap_to_tiles, LayersControl
from jupyter_functions import create_map
import warnings
warnings.filterwarnings("ignore")
# Setup display items
boundary_shp = os.path.join(output_folder,'geo_em.d01_boundary.shp')
b_shp = geopandas.read_file(boundary_shp)
b_shp = b_shp.to_crs(epsg=4326)
# Export vector to GeoJSON
b_json = os.path.join(output_folder, 'boundary.json')
b_shp.to_file(b_json, driver='GeoJSON')
# Read GeoJSON
with open(b_json, 'r') as f:
data = json.load(f)
# Obtain vector center point
x = b_shp.geometry.centroid.x
y = b_shp.geometry.centroid.y
map_center = y[0], x[0]
# Instantiate map object
m = Map(center=(41.50, -73.73), zoom=10, scroll_wheel_zoom=True)
# Read GeoJSON
with open(b_json, 'r') as f:
data = json.load(f)
# Obtain vector center point
x = b_shp.geometry.centroid.x
y = b_shp.geometry.centroid.y
map_center = y[0], x[0]
# Instantiate map object
m = create_map(map_center, 10)
# Read GeoJSON
geo_json = GeoJSON(data=data, name='Domain boundary')
# Define basemaps to swipe between
right_layer = basemap_to_tiles(basemap=basemaps.OpenStreetMap.Mapnik)
left_layer = basemap_to_tiles(basemap=basemaps.Esri.WorldImagery)
# Setup basemap swipe control
control = SplitMapControl(left_layer=left_layer, right_layer=right_layer)
m.add_control(control)
m.add_layer(geo_json)
# Draw map
m
The tool Build_GeoTiff_From_Geogrid_File.py
is a program to export variables from a WRF-Hydro input file (geogrid or Fulldom_hires) file to an output raster format, with all spatial and coordinate system metadata. If a 3-dimensional variable is selected, individual raster bands will be created in the output raster for each index in the 3rd dimension. If a 4-dimensional variable is selected, the first index in the 4th dimension will be selected and the variable will be treated as a 3-dimensional variable described above.
This tool is handy for performing a quick vizualization using GIS or othe software to examine the contents of the WRF-Hydro input file and overlay these grids with other goespatial data.
The tool takes three input parameters: an input Geogrid or Fulldom_hires netCDF file (-i
), a variable name (-v
), and an output GeoTiff raster file (-o
) that the tool will create. For this example, we will export the variable "HGT_M", or surface elevation in meters above sea level.
# Get script help information
! python Build_GeoTiff_From_Geogrid_File.py -h
Script initiated at Thu Nov 5 18:03:07 2020 usage: Build_GeoTiff_From_Geogrid_File.py [-h] -i IN_NC [-v VARIABLE] [-o OUT_FILE] This is a program to export >=2D variables from a WRF-Hydro input file (geogrid or Fulldom_hires) file to an output raster format, with all spatial and coordinate system metadata. If a 3-dimensional variable is selected, individual raster bands will be created in the output raster for each index in the 3rd dimension. If a 4-dimensional variable is selected, the first index in the 4th dimension will be selected and the variable will be treated as a 3-dimensional variable described above. optional arguments: -h, --help show this help message and exit -i IN_NC Path to WPS geogrid (geo_em.d0*.nc) file or WRF-Hydro Fulldom_hires.nc file. -v VARIABLE Name of the variable in the input netCDF file. default=HGT_M -o OUT_FILE Output GeoTiff raster file.
# Define the variable to export to raster
in_var = "HGT_M"
# Define the output raster file using variable name defined above
out_file = os.path.join(output_folder, f'{in_var}.tif')
# Print information to screen for reference
print('Command to run:\n')
print('python Build_GeoTiff_From_Geogrid_File.py \\\n\t -i {0} \\\n\t -v {1} \\\n\t -o {2}\n'.format(in_geogrid, in_var, out_file))
# Run the script with required parameters
! python Build_GeoTiff_From_Geogrid_File.py -i {in_geogrid} -v {in_var} -o {out_file}
Command to run: python Build_GeoTiff_From_Geogrid_File.py \ -i /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/geo_em.d01.nc \ -v HGT_M \ -o /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/HGT_M.tif Script initiated at Thu Nov 5 18:03:15 2020 Using default variable name: HGT_M Input WPS Geogrid or Fulldom file: /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/geo_em.d01.nc Input netCDF variable name: HGT_M Output raster file: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/HGT_M.tif WPS netCDF projection identification initiated... Map Projection: Lambert Conformal Conic Using MOAD_CEN_LAT for latitude of origin. Using Standard Parallel 2 in Lambert Conformal Conic map projection. Geo-referencing step completed without error in 0.05 seconds. X-dimension: 'west_east'. Y-dimension: 'south_north'. Reversing order of dimension 'south_north' Time dimension found: 'Time'. Time dimension size = 1. Dimensions and indices or slices on those dimensions: Time: 0 south_north: slice(None, None, -1) west_east: slice(None, None, None) Size of array being sent to raster: (16, 15) GDAL Data type derived from input array: 6 (float32) Bands in output raster: 1 Created GTiff format raster from HGT_M variable: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/HGT_M.tif Process complted in 0.19 seconds.
Now that the tool has completed, we want to take a look at the output. We will create another interactive map and load the data as a layer.
# Import third-party visualization libraries
import rasterio
from matplotlib import pyplot
from osgeo import gdal
from ipyleaflet import ImageOverlay
from jupyter_functions import cmap_options, show_raster_map
# Create a map object from pre-build function
m2 = create_map(map_center, 10)
# Render the map
m2
# Use pre-built function to render the GeoTiff on the map, already warped to the map's coordinate system
show_raster_map(out_file, m2, b_shp, output_folder)
Above, you will see the elevation grid applied to the map. This grid is 1km, so there is not much detail, but it is still useful to see if the topographic features are in the correct geographic locations according to the basemap. Also, there is no color-ramp for reference. This is a limitation of using the web browser over a GIS application.
The Build_Routing_Stack.py
script is a program to build the full set of hydrologically-processed routing grids and additional data required by WRF-Hydro. This is the main utility for performing WRF-Hydro GIS pre-processing. The required inputs are the domain file (WPS geogrid file), desired routing grid resolution as a function of geogrid resolution, and other options and parameter values. The output will be a "routing stack" zip file with WRF-Hydro domain and parameter files.
• Required Parameters:
-i
-WRF/WPS GEOGRID file (geo_em.d0*.nc)
-d
-High-resolution Elevation raster file (Esri GRID, GeoTIFF, VRT, etc.)
-R
-Regridding Factor – nesting relationship of routing:land surface model grid cells
-t
-Minimum basin area threshold (in routing grid cells)
-o
-Output ZIP File containing all script outputs
• Optional Parameters:
--CSV
-Station Locations location file (.csv)
-b
-Option to mask channel grids not contributing to provided station locations
-r
-Reach based (Muskingum / Muskingum-Cunge) routing option
-l
-Lake Polygons (polygon feature class or .shp)
-O
-OVROUGHRTFAC – Multiplier on Manning's roughness for overland flow. default=1.0
-T
-RETDEPRTFAC – Multiplier on maximum retention depth before flow is routed as overland flow. default=1.0
-LKSATFAC – (script global variable) Multiplier on saturated hydraulic conductivity in lateral flow direction. default=1000.0
--starts
-Path to point shapefile or feature class containing channel initiation locations (overrides -t
parameter)
--gw
-Path to polygon shapefile or feature class containing prescribed groundwater basin locations
This tool has many different parameters and possible configurations. Using the command line, we can take a look at the different arguements that can be used, if they are required or optional, and what their default values are.
! python Build_Routing_Stack.py -h
Script initiated at Thu Nov 5 18:04:19 2020 usage: Build_Routing_Stack.py [-h] -i IN_GEOGRID [--CSV IN_CSV] [-b BASIN_MASK] [-r RB_ROUTING] [-l IN_RESERVOIRS] -d INDEM [-R CELLSIZE] [-t THRESHOLD] [-o OUT_ZIP_FILE] [-O OVROUGHRTFAC_VAL] [-T RETDEPRTFAC_VAL] [--starts CHANNEL_STARTS] [--gw GW_POLYS] This is a program to perform the full routing-stack GIS pre-processingfor WRF- Hydro. The inputs will be related to the domain, the desired routing nest factor, and other options and parameter values. The output will be a routing stack zip file with WRF-Hydro domain and parameter files. optional arguments: -h, --help show this help message and exit -i IN_GEOGRID Path to WPS geogrid (geo_em.d0*.nc) file [REQUIRED] --CSV IN_CSV Path to input forecast point CSV file [OPTIONAL] -b BASIN_MASK Mask CHANNELGRID variable to forecast basins? [True/False]. default=False -r RB_ROUTING Create reach-based routing (RouteLink) files? [True/False]. default=False -l IN_RESERVOIRS Path to reservoirs shapefile or feature class [OPTIONAL]. If -l is TRUE, this is required. -d INDEM Path to input high-resolution elevation raster [REQUIRED] -R CELLSIZE Regridding (nest) Factor. default=10 -t THRESHOLD Number of routing grid cells to define stream. default=200 -o OUT_ZIP_FILE Output routing stack ZIP file -O OVROUGHRTFAC_VAL OVROUGHRTFAC value. default=1.0 -T RETDEPRTFAC_VAL RETDEPRTFAC value. default=1.0 --starts CHANNEL_STARTS Path to channel initiation points feature class. Must be 2D point type. [OPTIONAL] --gw GW_POLYS Path to groundwater polygons feature class [OPTIONAL]
We will begin by assigning file paths to python variables, then substitude those values in the command line syntax.
import Build_Routing_Stack
# Define script input parameters using python variables
in_geogrid = os.path.join(data_folder, 'geo_em.d01.nc')
lakes = os.path.join(data_folder, 'lake_shapes', 'lakes.shp')
csv = os.path.join(data_folder, 'croton_frxst_pts_FOSS.csv')
in_dem = os.path.join(data_folder, 'NED_30m_Croton.tif')
regrid_factor = 4
routing_cells = 20
out_zip = os.path.join(output_folder, 'croton_test.zip')
# Print information to screen for reference
print('Command to run:\n')
print('python Build_Routing_Stack.py \\\n\t -i {0} \\\n\t -l {1} \\\n\t --CSV {2} \\\n\t -d {3} \\\n\t -R {4} \\\n\t -t {5} \\\n\t -o {6}\n'.format(in_geogrid, lakes, csv, in_dem, regrid_factor, routing_cells, out_zip))
# Run the script with required parameters
! python Build_Routing_Stack.py -i {in_geogrid} -l {lakes} --CSV {csv} -d {in_dem} -R {regrid_factor} -t {routing_cells} -o {out_zip}
Command to run: python Build_Routing_Stack.py \ -i /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/geo_em.d01.nc \ -l /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/lake_shapes/lakes.shp \ --CSV /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv \ -d /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/NED_30m_Croton.tif \ -R 4 \ -t 20 \ -o /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/croton_test.zip Script initiated at Thu Nov 5 18:04:26 2020 Parameter values that have not been altered from script default values: Using default basin mask setting: False Using default reach-based routing setting: False Using default OVROUGHRTFAC parameter value: 1.0 Using default RETDEPRTFAC parameter value: 1.0 Values that will be used in building this routing stack: Input WPS Geogrid file: /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/geo_em.d01.nc Forecast Point CSV file: /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv Mask CHANNELGRID variable to forecast basins?: False Create reach-based routing (RouteLink) files?: False Lake polygon feature class: /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/lake_shapes/lakes.shp Input high-resolution DEM: /home/docker/.mnt/common/mollymca/GIS_Training/Croton_Lambert/NED_30m_Croton.tif Regridding factor: 4 Stream initiation threshold: 20 OVROUGHRTFAC parameter value: 1.0 RETDEPRTFAC parameter value: 1.0 Input channel initiation start point feature class: None Input groundwater basin polygons: None Output ZIP file: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/croton_test.zip Running Process GEOGRID function Forecast points provided. Reach-based routing files will not be created. WPS netCDF projection identification initiated... Map Projection: Lambert Conformal Conic Using MOAD_CEN_LAT for latitude of origin. Using Standard Parallel 2 in Lambert Conformal Conic map projection. Geo-referencing step completed without error in 0.05 seconds. Building sub-grid of model grid. Original grid spacing dx=1000.0, dy=-1000.0 Original grid size: rows=16, cols=15 New grid spacing: dx=250.0, dy=-250.0 New dimensions: rows=64, cols=60 Created projection definition from input NetCDF GEOGRID file. Proj4: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs Coarse grid GeoTransform: 1841999.0194101864 1000.0 0 278495.814356732 0 -1000.0 Coarse grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732] Fine grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732] GDAL Data type derived from input array: 6 (float32) Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. netCDF global attributes set after 0.07 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.31 seconds. Deriving geocentric coordinates on routing grid from bilinear interpolation of geogrid coordinates. GDAL Data type derived from input array: 6 (float32) GDAL Data type derived from input array: 6 (float32) Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. Proceeding to add LATITUDE and LONGITUDE variables after 0.14 seconds. netCDF global attributes set after 0.14 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Process: landuse written to output netCDF. Terrain processing step initiated... Using WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020 Depression Filling algorithm: Whitebox Fill Depressions. Process: TOPOGRAPHY written to output netCDF. Coerced 21 0-value flow direction cells to flow off of the grid. Process: FLOWDIRECTION written to output netCDF. Process: FLOWACC written to output netCDF. Flow accumulation will be thresholded to build channel pixels. Process: CHANNELGRID written to output netCDF. Process: STREAMORDER written to output netCDF. Process: RETDEPRTFAC written to output netCDF. Process: OVROUGHRTFAC written to output netCDF. Process: LKSATFAC written to output netCDF. Terrain processing step completed without error in 0.28 seconds. Forecast points provided and basins being delineated. Process: frxst_pts written to output netCDF. Process: basn_msk written to output netCDF. Channelgrid will not be masked to basins. Built forecast point outputs in 0.45 seconds. Reservoir polygons provided. Lake routing will be activated. Adding reservoirs to routing stack. Gridded: True MEMORY driver is available. Done producing output vector polygon shapefile in 0.00 seconds Input shapefile projection does not match requested output. Transforming. Number of output features: 1 of 1 Completed reprojection and-or clipping in 0.17 seconds. Adding auto-incremented lake ID field (1...n) Starting to gather lake centroid and area information. Done gathering lake centroid information. Found 1 lakes on active channels. Lost 0 lakes that were not on active channels. Removing lakes not on gridded channel network Process: LAKEGRID written to output netCDF. Process: CHANNELGRID written to output netCDF. Starting to create lake parameter table. Lakes Table: 1 Lakes Starting to fill in lake parameter table NC file. Done writing LAKEPARM.nc table to disk. Lake parameter table created without error in 0.78 seconds. Beginning to build 2D groundwater basin inputs Building groundwater inputs using FullDom LINKID local basins Generating LINKID grid for building local sub-basins. Finished building fine-grid groundwater basin grids in 0.10 seconds Beginning to build coarse-grid groundwater basins and parameters Found 72 basins in the watershed grid Raster resampling initiated... The High-resolution dataset will be 1000.0m Projected input raster to model grid in 0.02 seconds. Found 69 basins (potentially including nodata values) in the file after resampling to the coarse grid. GDAL Data type derived from input array: 5 (int32) Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. netCDF global attributes set after 0.06 seconds. NC dimensions: 16, 15 GWBUCKS array dimensions: 16, 15 Process: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/scratchdir/GWBASINS.nc completed without error Finished building groundwater grid file in 0.11 seconds Calculating size and ID parameters for basin polygons. Created output bucket parameter table (.nc): /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/scratchdir/GWBUCKPARM.nc. Finished building groundwater bucket parameter table in 0.06 seconds. Finished building groundwater parameter files in 0.18 seconds Built output .zip file in 2.76 seconds. Process completed in 2.87 seconds.
The Build_Routing_Stack.py
script creates a Zip archive of output files according to the options provided to the tool. There will be at least four netCDF files. The output Zip file may additionally include shapefiles (.shp and accompanying files) describing the geometry of modeled lakes or the vector stream network. Below is an alphabetically sorted list of gridded variables that are created by the Build_Routing_Stack.py
tool.
Fulldom_hires.nc file. This file stores 2D gridded variables that describe the hydro routing grid:
-b
option is set to TRUE, the output will be masked to the gaged basins provided, where non-gaged channels are given a value of ‘-1’. If lake routing is activated, lake outflow points will be identified by the lake ID value.-l
parameter, this grid will contain ID values for each lake that can be resolved on the routing grid. Otherwise, this grid will be uniform with values of -9999.-r
is TRUE.-O
parameter.-T
parameter.-d INDEM
), which should be specified in meters (m) above sea level (ASL). This grid is derived from the elevation values in the input elevation raster, but has been resampled to the routing grid resolution, and pit-filled to remove depressions.Other files:
-l
parameter. The table will contain a record for each lake feature in the + Fulldom_hires.nc ‘LAKEGRID’ variable, and contain derived and default parameters for each lake.-r
parameter is TRUE. The file contains a record for each stream segment. The stream segments in this table are also identified by the unique ‘LINKID’ values in the ‘LINKID’ variable in the ‘Fulldom_hires.nc’ file, and values in the ‘STRM_VAL’ field of the output ‘streams.shp’ shapefile. This table contains derived and default stream segment parameters that are calculated based on the vector stream network and topology in the ‘streams.shp’ shapefile.-l
is used. The ‘STRM_VAL’ field is the unique identifier for each stream segment, and corresponds to the ‘link’ variable of the ‘Route_Link.nc’ file and the ‘LINKID’ variable in the ‘Fulldom_hires.nc’ file. The geometry of the stream segments in this shapefile informs many of the parameters in the ‘Route_Link.nc’ file.-r TRUE
is used, then Fulldom_hires.nc 'LAKEID' variable will contain -9999 values only. The geometry of reservoir objects informs many of the parameters in 'LAKEPARM.nc' file.The Examine_Outputs_of_GIS_Preprocessor.py
script is a tool used to create easy-to-use geospatial files for examining the resulting WRF-Hydro input files in a GIS. This script takes the output ZIP file generated using the Process Geogrid script (executed above) and creates a raster from each 2D variable in each WRF-Hydro input netCDF file. In addition, other data will be extracted such as any shapefiles, 1D netCDF tables, etc. The input to the tool should be a .zip file that was created using the WRF Hydro pre-processing tools. The tool will create the output folder if it does not already exist, and write all results to that location.
This tool has a single input and single output parameter. Using the command line, we can take a look at the arguments.
! python Examine_Outputs_of_GIS_Preprocessor.py -h
Script initiated at Thu Nov 5 18:04:34 2020 usage: Examine_Outputs_of_GIS_Preprocessor.py [-h] -i IN_ZIP -o OUT_FOLDER This tool takes the output zip file from the ProcessGeogrid script and creates a raster from each output NetCDF file. The Input should be a .zip file that was created using theWRF Hydro pre-processing tools. The tool will create the folder which will contain theresults (out_folder), if that folder does not already exist. optional arguments: -h, --help show this help message and exit -i IN_ZIP Path to WRF Hydro routing grids zip file. -o OUT_FOLDER Path to output folder.
We will define the output directory for the tool to create, then use defined variables to execute the command line tool.
# Define output directory to store GeoTiff output of all routing stack grids
raster_outputs = os.path.join(output_folder, "Raster_Outputs")
# Print information to screen for reference
print('Command to run:\n')
print('python Examine_Outputs_of_GIS_Preprocessor.py \\\n\t -i {0} \\\n\t -o {1}\n'.format(out_zip, raster_outputs))
# Run the script with required parameters
! python Examine_Outputs_of_GIS_Preprocessor.py -i {out_zip} -o {raster_outputs}
Command to run: python Examine_Outputs_of_GIS_Preprocessor.py \ -i /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/croton_test.zip \ -o /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs Script initiated at Thu Nov 5 18:04:41 2020 File Copied: GEOGRID_LDASOUT_Spatial_Metadata.nc File Copied: LAKEPARM.nc GeoTransform: (1841999.0194101864, 250.0, 0.0, 278495.814356732, 0.0, -250.0) DX: 250.0 DY: 250.0 PROJ.4 string: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/CHANNELGRID.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/FLOWDIRECTION.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/FLOWACC.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/TOPOGRAPHY.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/RETDEPRTFAC.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/OVROUGHRTFAC.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/STREAMORDER.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/frxst_pts.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/basn_msk.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/LAKEGRID.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/landuse.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/LKSATFAC.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/LATITUDE.tif File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/LONGITUDE.tif File Copied: GWBUCKPARM.nc File Copied: lakes.dbf File Copied: lakes.prj File Copied: lakes.shx File Copied: lakes.shp GeoTransform: (1841999.0194101864, 1000.0, 0.0, 278495.814356732, 0.0, -1000.0) DX: 1000.0 DY: 1000.0 PROJ.4 string: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs File Created: /home/docker/.mnt/common/mollymca/GIS_Training/Outputs/Raster_Outputs/BASIN.tif Extraction of WRF routing grids completed. Process complted in 1.81 seconds.
Finally, we want to take a look at some of the output rasters. Below, we utilize a Jupyter widget to choose each raster from a drop down menu. Each of the 2D variables in Fulldom_hires.nc (the routing grid netCDF file) will be displayed as a rectangular grid. Take a look at some of the outputs.
import ipywidgets as widgets
from ipywidgets import interact
def see_raster(x):
src = rasterio.open(os.path.join(raster_outputs, f"{x}.tif"))
cmap, norm = cmap_options(x)
if x in ['TOPOGRAPHY']:
pyplot.imshow(src.read(1), cmap=cmap, aspect='auto', norm=norm, interpolation='nearest', vmin=0)
else:
pyplot.imshow(src.read(1), cmap=cmap, aspect='auto', norm=norm, interpolation='nearest')
cbar = pyplot.colorbar()
# Keep the automatic aspect while scaling the image up in size
fig = pyplot.gcf()
w, h = fig.get_size_inches()
fig.set_size_inches(w * 1.75, h * 1.75)
# Show image
pyplot.show()
in_raster = widgets.Dropdown(
options=[('Basin', 'BASIN'), ('Basin mask', 'basn_msk'), ('Channel grid', 'CHANNELGRID'), ('Flow accumulation', 'FLOWACC'),
('Flow direction', 'FLOWDIRECTION'), ('Forecast points', 'frxst_pts'), ('Lake grid', 'LAKEGRID'),
('Land use', 'landuse'), ('Latitude', 'LATITUDE'), ('LKSATFAC', 'LKSATFAC'), ('Longitude', 'LONGITUDE'),
('OVROUGHRTFAC', 'OVROUGHRTFAC'), ('RETDEPRTFAC', 'RETDEPRTFAC'), ('Stream order', 'STREAMORDER'),
('Topography', 'TOPOGRAPHY')],
value='FLOWACC',
description='Variable:')
interact(see_raster, x=in_raster)
<function __main__.see_raster(x)>
We can also look at the raster data projected onto a map. Use the dropdown widget below the map to add rasters as layers to the map. Use the menu in the top right corner of the map to turn layers on and off. Behind the scenes, each model grid is being reprojected to Web Mercator for overlay on a WMS basemap for some interactivity.
# Create a map object from pre-build function
m3 = create_map(map_center, 10) # Create the map
m3 # Render the map
def f(x):
r_path = os.path.join(raster_outputs, f"{x}.tif")
show_raster_map(r_path, m3, b_shp, output_folder)
in_raster = widgets.Dropdown(
options=[('Basin', 'BASIN'), ('Forecast point basins', 'basn_msk'), ('Channel grid', 'CHANNELGRID'), ('Flow accumulation', 'FLOWACC'), ('Flow direction', 'FLOWDIRECTION'),
('Forecast points', 'frxst_pts'), ('Lake grid', 'LAKEGRID'), ('Landuse', 'landuse'), ('OVROUGHRTFAC grid', 'OVROUGHRTFAC'), ('RETDEPRTFAC grid', 'RETDEPRTFAC'),
('Stream order', 'STREAMORDER'), ('Topography', 'TOPOGRAPHY')],
value='TOPOGRAPHY',
description='Variable:',
)
interact(f, x=in_raster)
<function __main__.f(x)>
We will run command-line arguments to create output directories for each of the test cases. Then, each cell will run the WRF-Hydro GIS Pre-processor to build the 'routing stack' files. The final command is used to unzip the outputs so they can be more readily used by WRF-Hydro. There are 4 Non-NWM test-case configurations which can all be created using different combinations of arguments to the GIS pre-processor. Currently missing at this time is the Reach-with-Lakes capability as well as the ability to build wrfinput, hydro2dtbl, and soil_properties files.
This domain is a gridded domain with lakes, gages and a regridding factor of 4 and threshold of 20
# Create the output directory and ensure it is empty
! mkdir -p /home/docker/GIS_Training/Outputs/Gridded
! rm -rf /home/docker/GIS_Training/Outputs/Gridded/*
# Run the GIS Pre-processing script (with line-breaks)
! python Build_Routing_Stack.py \
-i /home/docker/wrf-hydro-training/example_case/Gridded/DOMAIN/geo_em.d01.nc \
-l /home/docker/GIS_Training/Croton_Lambert/lake_shapes/lakes.shp \
--CSV /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv \
-d /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif \
-R 4 \
-t 20 \
-o /home/docker/GIS_Training/Outputs/Gridded/Gridded_r4_t20_lakes_frxst_mask.zip
# Unzip the directory in-place
! unzip /home/docker/GIS_Training/Outputs/Gridded/Gridded_r4_t20_lakes_frxst_mask.zip -d /home/docker/GIS_Training/Outputs/Gridded
Script initiated at Thu Nov 5 18:05:35 2020 Parameter values that have not been altered from script default values: Using default basin mask setting: False Using default reach-based routing setting: False Using default OVROUGHRTFAC parameter value: 1.0 Using default RETDEPRTFAC parameter value: 1.0 Values that will be used in building this routing stack: Input WPS Geogrid file: /home/docker/wrf-hydro-training/example_case/Gridded/DOMAIN/geo_em.d01.nc Forecast Point CSV file: /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv Mask CHANNELGRID variable to forecast basins?: False Create reach-based routing (RouteLink) files?: False Lake polygon feature class: /home/docker/GIS_Training/Croton_Lambert/lake_shapes/lakes.shp Input high-resolution DEM: /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif Regridding factor: 4 Stream initiation threshold: 20 OVROUGHRTFAC parameter value: 1.0 RETDEPRTFAC parameter value: 1.0 Input channel initiation start point feature class: None Input groundwater basin polygons: None Output ZIP file: /home/docker/GIS_Training/Outputs/Gridded/Gridded_r4_t20_lakes_frxst_mask.zip Running Process GEOGRID function Forecast points provided. Reach-based routing files will not be created. WPS netCDF projection identification initiated... Map Projection: Lambert Conformal Conic Using MOAD_CEN_LAT for latitude of origin. Using Standard Parallel 2 in Lambert Conformal Conic map projection. Geo-referencing step completed without error in 0.05 seconds. Building sub-grid of model grid. Original grid spacing dx=1000.0, dy=-1000.0 Original grid size: rows=16, cols=15 New grid spacing: dx=250.0, dy=-250.0 New dimensions: rows=64, cols=60 Created projection definition from input NetCDF GEOGRID file. Proj4: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs Coarse grid GeoTransform: 1841999.0194101864 1000.0 0 278495.814356732 0 -1000.0 Coarse grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732] Fine grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732] GDAL Data type derived from input array: 6 (float32) Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. netCDF global attributes set after 0.14 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.23 seconds. Deriving geocentric coordinates on routing grid from bilinear interpolation of geogrid coordinates. GDAL Data type derived from input array: 6 (float32) GDAL Data type derived from input array: 6 (float32) Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. Proceeding to add LATITUDE and LONGITUDE variables after 0.19 seconds. netCDF global attributes set after 0.20 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Process: landuse written to output netCDF. Terrain processing step initiated... Using WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020 Depression Filling algorithm: Whitebox Fill Depressions. Process: TOPOGRAPHY written to output netCDF. Coerced 21 0-value flow direction cells to flow off of the grid. Process: FLOWDIRECTION written to output netCDF. Process: FLOWACC written to output netCDF. Flow accumulation will be thresholded to build channel pixels. Process: CHANNELGRID written to output netCDF. Process: STREAMORDER written to output netCDF. Process: RETDEPRTFAC written to output netCDF. Process: OVROUGHRTFAC written to output netCDF. Process: LKSATFAC written to output netCDF. Terrain processing step completed without error in 0.27 seconds. Forecast points provided and basins being delineated. Process: frxst_pts written to output netCDF. Process: basn_msk written to output netCDF. Channelgrid will not be masked to basins. Built forecast point outputs in 0.55 seconds. Reservoir polygons provided. Lake routing will be activated. Adding reservoirs to routing stack. Gridded: True MEMORY driver is available. Done producing output vector polygon shapefile in 0.00 seconds Input shapefile projection does not match requested output. Transforming. Number of output features: 1 of 1 Completed reprojection and-or clipping in 0.16 seconds. Adding auto-incremented lake ID field (1...n) Starting to gather lake centroid and area information. Done gathering lake centroid information. Found 1 lakes on active channels. Lost 0 lakes that were not on active channels. Removing lakes not on gridded channel network Process: LAKEGRID written to output netCDF. Process: CHANNELGRID written to output netCDF. Starting to create lake parameter table. Lakes Table: 1 Lakes Starting to fill in lake parameter table NC file. Done writing LAKEPARM.nc table to disk. Lake parameter table created without error in 0.94 seconds. Beginning to build 2D groundwater basin inputs Building groundwater inputs using FullDom LINKID local basins Generating LINKID grid for building local sub-basins. Finished building fine-grid groundwater basin grids in 0.28 seconds Beginning to build coarse-grid groundwater basins and parameters Found 72 basins in the watershed grid Raster resampling initiated... The High-resolution dataset will be 1000.0m Projected input raster to model grid in 0.02 seconds. Found 69 basins (potentially including nodata values) in the file after resampling to the coarse grid. GDAL Data type derived from input array: 5 (int32) Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. netCDF global attributes set after 0.07 seconds. NC dimensions: 16, 15 GWBUCKS array dimensions: 16, 15 Process: /home/docker/GIS_Training/Outputs/Gridded/scratchdir/GWBASINS.nc completed without error Finished building groundwater grid file in 0.24 seconds Calculating size and ID parameters for basin polygons. Created output bucket parameter table (.nc): /home/docker/GIS_Training/Outputs/Gridded/scratchdir/GWBUCKPARM.nc. Finished building groundwater bucket parameter table in 0.11 seconds. Finished building groundwater parameter files in 0.36 seconds Built output .zip file in 3.52 seconds. Process completed in 3.68 seconds. Archive: /home/docker/GIS_Training/Outputs/Gridded/Gridded_r4_t20_lakes_frxst_mask.zip inflating: /home/docker/GIS_Training/Outputs/Gridded/GEOGRID_LDASOUT_Spatial_Metadata.nc inflating: /home/docker/GIS_Training/Outputs/Gridded/LAKEPARM.nc inflating: /home/docker/GIS_Training/Outputs/Gridded/Fulldom_hires.nc inflating: /home/docker/GIS_Training/Outputs/Gridded/GWBUCKPARM.nc inflating: /home/docker/GIS_Training/Outputs/Gridded/lakes.dbf inflating: /home/docker/GIS_Training/Outputs/Gridded/lakes.prj inflating: /home/docker/GIS_Training/Outputs/Gridded/lakes.shx inflating: /home/docker/GIS_Training/Outputs/Gridded/lakes.shp inflating: /home/docker/GIS_Training/Outputs/Gridded/GWBASINS.nc
This domain is a gridded domain without lakes. It contains gages, masked basins and a regridding factor of 4 and thrshold of 20
# Create the output directory and ensure it is empty
! mkdir -p /home/docker/GIS_Training/Outputs/Gridded_no_lakes
! rm -rf /home/docker/GIS_Training/Outputs/Gridded_no_lakes/*
# Run the GIS Pre-processing script (with line-breaks)
! python Build_Routing_Stack.py \
-i /home/docker/wrf-hydro-training/example_case/Gridded_no_lakes/DOMAIN/geo_em.d01.nc \
--CSV /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv \
-d /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif \
-R 4 \
-t 20 \
-o /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Gridded_r4_t20_frxst_mask.zip
# Unzip the directory in-place
! unzip /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Gridded_r4_t20_frxst_mask.zip -d /home/docker/GIS_Training/Outputs/Gridded_no_lakes/
Script initiated at Thu Nov 5 18:06:56 2020 Parameter values that have not been altered from script default values: Using default basin mask setting: False Using default reach-based routing setting: False Using default OVROUGHRTFAC parameter value: 1.0 Using default RETDEPRTFAC parameter value: 1.0 Values that will be used in building this routing stack: Input WPS Geogrid file: /home/docker/wrf-hydro-training/example_case/Gridded_no_lakes/DOMAIN/geo_em.d01.nc Forecast Point CSV file: /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv Mask CHANNELGRID variable to forecast basins?: False Create reach-based routing (RouteLink) files?: False Lake polygon feature class: Input high-resolution DEM: /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif Regridding factor: 4 Stream initiation threshold: 20 OVROUGHRTFAC parameter value: 1.0 RETDEPRTFAC parameter value: 1.0 Input channel initiation start point feature class: None Input groundwater basin polygons: None Output ZIP file: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Gridded_r4_t20_frxst_mask.zip Running Process GEOGRID function Forecast points provided. Reach-based routing files will not be created. WPS netCDF projection identification initiated... Map Projection: Lambert Conformal Conic Using MOAD_CEN_LAT for latitude of origin. Using Standard Parallel 2 in Lambert Conformal Conic map projection. Geo-referencing step completed without error in 0.05 seconds. Building sub-grid of model grid. Original grid spacing dx=1000.0, dy=-1000.0 Original grid size: rows=16, cols=15 New grid spacing: dx=250.0, dy=-250.0 New dimensions: rows=64, cols=60 Created projection definition from input NetCDF GEOGRID file. Proj4: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs Coarse grid GeoTransform: 1841999.0194101864 1000.0 0 278495.814356732 0 -1000.0 Coarse grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732] Fine grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732] GDAL Data type derived from input array: 6 (float32) Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. netCDF global attributes set after 0.06 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.27 seconds. Deriving geocentric coordinates on routing grid from bilinear interpolation of geogrid coordinates. GDAL Data type derived from input array: 6 (float32) GDAL Data type derived from input array: 6 (float32) Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. Proceeding to add LATITUDE and LONGITUDE variables after 0.47 seconds. netCDF global attributes set after 0.47 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Process: landuse written to output netCDF. Terrain processing step initiated... Using WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020 Depression Filling algorithm: Whitebox Fill Depressions. Process: TOPOGRAPHY written to output netCDF. Coerced 21 0-value flow direction cells to flow off of the grid. Process: FLOWDIRECTION written to output netCDF. Process: FLOWACC written to output netCDF. Flow accumulation will be thresholded to build channel pixels. Process: CHANNELGRID written to output netCDF. Process: STREAMORDER written to output netCDF. Process: RETDEPRTFAC written to output netCDF. Process: OVROUGHRTFAC written to output netCDF. Process: LKSATFAC written to output netCDF. Terrain processing step completed without error in 0.27 seconds. Forecast points provided and basins being delineated. Process: frxst_pts written to output netCDF. Process: basn_msk written to output netCDF. Channelgrid will not be masked to basins. Built forecast point outputs in 0.51 seconds. Beginning to build 2D groundwater basin inputs Building groundwater inputs using FullDom LINKID local basins Generating LINKID grid for building local sub-basins. Finished building fine-grid groundwater basin grids in 0.11 seconds Beginning to build coarse-grid groundwater basins and parameters Found 72 basins in the watershed grid Raster resampling initiated... The High-resolution dataset will be 1000.0m Projected input raster to model grid in 0.02 seconds. Found 69 basins (potentially including nodata values) in the file after resampling to the coarse grid. GDAL Data type derived from input array: 5 (int32) Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. netCDF global attributes set after 0.06 seconds. NC dimensions: 16, 15 GWBUCKS array dimensions: 16, 15 Process: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/scratchdir/GWBASINS.nc completed without error Finished building groundwater grid file in 0.15 seconds Calculating size and ID parameters for basin polygons. Created output bucket parameter table (.nc): /home/docker/GIS_Training/Outputs/Gridded_no_lakes/scratchdir/GWBUCKPARM.nc. Finished building groundwater bucket parameter table in 0.08 seconds. Finished building groundwater parameter files in 0.25 seconds Built output .zip file in 2.46 seconds. Process completed in 2.57 seconds. Archive: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Gridded_r4_t20_frxst_mask.zip inflating: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/GEOGRID_LDASOUT_Spatial_Metadata.nc inflating: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/Fulldom_hires.nc inflating: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/GWBUCKPARM.nc inflating: /home/docker/GIS_Training/Outputs/Gridded_no_lakes/GWBASINS.nc
This domain is a reach-based routing configuration with gages, no masking, no lakes. Regridding factor of 4 and threshold of 20
# Create the output directory and ensure it is empty
! mkdir -p /home/docker/GIS_Training/Outputs/Reach
! rm -rf /home/docker/GIS_Training/Outputs/Reach/*
# Run the GIS Pre-processing script (with line-breaks)
! python Build_Routing_Stack.py \
-i /home/docker/wrf-hydro-training/example_case/Reach/DOMAIN/geo_em.d01.nc \
-r True \
--CSV /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv \
-d /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif \
-R 4 \
-t 20 \
-o /home/docker/GIS_Training/Outputs/Reach/Reach_r4_t20_frxst.zip
# Unzip the directory in-place
! unzip /home/docker/GIS_Training/Outputs/Reach/Reach_r4_t20_frxst.zip -d /home/docker/GIS_Training/Outputs/Reach/
Script initiated at Thu Nov 5 18:07:03 2020 Parameter values that have not been altered from script default values: Using default basin mask setting: False Using default OVROUGHRTFAC parameter value: 1.0 Using default RETDEPRTFAC parameter value: 1.0 Values that will be used in building this routing stack: Input WPS Geogrid file: /home/docker/wrf-hydro-training/example_case/Reach/DOMAIN/geo_em.d01.nc Forecast Point CSV file: /home/docker/GIS_Training/Croton_Lambert/croton_frxst_pts_FOSS.csv Mask CHANNELGRID variable to forecast basins?: False Create reach-based routing (RouteLink) files?: True Lake polygon feature class: Input high-resolution DEM: /home/docker/GIS_Training/Croton_Lambert/NED_30m_Croton.tif Regridding factor: 4 Stream initiation threshold: 20 OVROUGHRTFAC parameter value: 1.0 RETDEPRTFAC parameter value: 1.0 Input channel initiation start point feature class: None Input groundwater basin polygons: None Output ZIP file: /home/docker/GIS_Training/Outputs/Reach/Reach_r4_t20_frxst.zip Running Process GEOGRID function Forecast points provided. Reach-based routing files will be created. WPS netCDF projection identification initiated... Map Projection: Lambert Conformal Conic Using MOAD_CEN_LAT for latitude of origin. Using Standard Parallel 2 in Lambert Conformal Conic map projection. Geo-referencing step completed without error in 0.05 seconds. Building sub-grid of model grid. Original grid spacing dx=1000.0, dy=-1000.0 Original grid size: rows=16, cols=15 New grid spacing: dx=250.0, dy=-250.0 New dimensions: rows=64, cols=60 Created projection definition from input NetCDF GEOGRID file. Proj4: +proj=lcc +lat_0=41.4710083007812 +lon_0=-97 +lat_1=30 +lat_2=60 +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs Coarse grid GeoTransform: 1841999.0194101864 1000.0 0 278495.814356732 0 -1000.0 Coarse grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732] Fine grid extent [Xmin, Ymin, Xmax, Ymax]: [1841999.0194101864, 262495.814356732, 1856999.0194101864, 278495.814356732] GDAL Data type derived from input array: 6 (float32) Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. netCDF global attributes set after 0.12 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.23 seconds. Deriving geocentric coordinates on routing grid from bilinear interpolation of geogrid coordinates. GDAL Data type derived from input array: 6 (float32) GDAL Data type derived from input array: 6 (float32) Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. Proceeding to add LATITUDE and LONGITUDE variables after 0.15 seconds. netCDF global attributes set after 0.15 seconds. Raster resampling initiated... The High-resolution dataset will be 250.0m Projected input raster to model grid in 0.02 seconds. Process: landuse written to output netCDF. Terrain processing step initiated... Using WhiteboxTools v1.4.0 by Dr. John B. Lindsay (c) 2017-2020 Depression Filling algorithm: Whitebox Fill Depressions. Process: TOPOGRAPHY written to output netCDF. Coerced 21 0-value flow direction cells to flow off of the grid. Process: FLOWDIRECTION written to output netCDF. Process: FLOWACC written to output netCDF. Flow accumulation will be thresholded to build channel pixels. Process: CHANNELGRID written to output netCDF. Process: STREAMORDER written to output netCDF. Process: RETDEPRTFAC written to output netCDF. Process: OVROUGHRTFAC written to output netCDF. Process: LKSATFAC written to output netCDF. Terrain processing step completed without error in 0.36 seconds. Forecast points provided and basins being delineated. Process: frxst_pts written to output netCDF. Process: basn_msk written to output netCDF. Channelgrid will not be masked to basins. Built forecast point outputs in 0.55 seconds. Routing table will be created... Stream to features step complete. Warning: Number of unique IDs exceeds limit of 16-bit unsigned integer type. Not all reaches may be converted to stream vectors. Check output carefully. Found 69 unique IDs in stream vector layer. Eliminating 3 IDs in LINKID grid that could not be resolved in stream vector layer. [43, 44, 65] Adding forecast points:LINKID association. Found 3 forecast point:LINKID associations. All dictionaries have been created in 0.00 seconds. Fields have been added to the shapefile. Projection defined for input vector layer in 0.13 seconds. Starting to fill in routing table NC file. Done writing NC file to disk. Routing table created without error. Reach-based routing inputs generated in 0.78 seconds. Beginning to build 2D groundwater basin inputs Building groundwater inputs using FullDom LINKID local basins Generating LINKID grid for building local sub-basins. Finished building fine-grid groundwater basin grids in 0.07 seconds Beginning to build coarse-grid groundwater basins and parameters Found 72 basins in the watershed grid Raster resampling initiated... The High-resolution dataset will be 1000.0m Projected input raster to model grid in 0.02 seconds. Found 69 basins (potentially including nodata values) in the file after resampling to the coarse grid. GDAL Data type derived from input array: 5 (int32) Creating CF-netCDF File. Esri PE String: PROJCS["unnamed",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["unnamed",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-97.0],PARAMETER["Standard_Parallel_1",30.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Latitude_Of_Origin",41.4710083007812],UNIT["Meter",1.0]] Map Projection of input raster : lambert_conformal_conic Starting Process: Building to XMap/YMap Conversion of input raster to XMap/YMap completed without error. netCDF global attributes set after 0.14 seconds. NC dimensions: 16, 15 GWBUCKS array dimensions: 16, 15 Process: /home/docker/GIS_Training/Outputs/Reach/scratchdir/GWBASINS.nc completed without error Finished building groundwater grid file in 0.17 seconds Calculating size and ID parameters for basin polygons. Created output bucket parameter table (.nc): /home/docker/GIS_Training/Outputs/Reach/scratchdir/GWBUCKPARM.nc. Finished building groundwater bucket parameter table in 0.11 seconds. Finished building groundwater parameter files in 0.30 seconds Built output .zip file in 3.14 seconds. Process completed in 3.23 seconds. Archive: /home/docker/GIS_Training/Outputs/Reach/Reach_r4_t20_frxst.zip inflating: /home/docker/GIS_Training/Outputs/Reach/streams.shx inflating: /home/docker/GIS_Training/Outputs/Reach/streams.shp inflating: /home/docker/GIS_Training/Outputs/Reach/GEOGRID_LDASOUT_Spatial_Metadata.nc inflating: /home/docker/GIS_Training/Outputs/Reach/Fulldom_hires.nc inflating: /home/docker/GIS_Training/Outputs/Reach/GWBUCKPARM.nc inflating: /home/docker/GIS_Training/Outputs/Reach/streams.dbf inflating: /home/docker/GIS_Training/Outputs/Reach/streams.prj inflating: /home/docker/GIS_Training/Outputs/Reach/GWBASINS.nc inflating: /home/docker/GIS_Training/Outputs/Reach/Route_Link.nc