Skip to content
This repository has been archived by the owner on May 20, 2024. It is now read-only.

Commit

Permalink
subset and output bbox for CONUS1.0 domain
Browse files Browse the repository at this point in the history
  • Loading branch information
@hoangtv1899
hoangtv1899 committed Apr 8, 2020
1 parent c11f0e3 commit abcdf3c
Show file tree
Hide file tree
Showing 4 changed files with 120 additions and 377 deletions.
83 changes: 3 additions & 80 deletions Clip_Inputs/clip_inputs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,33 +125,10 @@ def subset(arr,mask_arr,ds_ref, crop_to_domain, ndata=0):
#parser_a.add_argument('-z_bottom',type=int, help = 'bottom of domain (optional). Default is 0')
#parser_a.add_argument('-z_top',type=int, help = 'top of domain (optional). Default is 1000')

#group 2: using mask file
parser_b = subparsers.add_parser('mask', help='subset using a mask file')
parser_b.add_argument('-mask_file',type=str, help = 'input mask file')
parser_b.add_argument('-out_name',type=str, help = 'name of output solidfile (optional)')
parser_b.add_argument('-dx',type=int, help = 'spatial resolution of solidfile (optional). Default is 1000')
parser_b.add_argument('-dz',type=int, help = 'lateral resolution of solidfile (optional). Default is 1000')
parser_b.add_argument('-printmask',type=int, help = 'print mask (optional). Default is 0')
parser_b.add_argument('-printbbox',type=int, help = 'print bounding box (optional). Default is 0')
#parser_b.add_argument('-z_bottom',type=int, help = 'bottom of domain (optional). Default is 0')
#parser_b.add_argument('-z_top',type=int, help = 'top of domain (optional). Default is 1000')

#group 3: using custom watershed
parser_c = subparsers.add_parser('define_watershed', help='subset using a newly created watershed')
parser_c.add_argument('-dir_file',type=str, help = 'input direction file',)
parser_c.add_argument('-outlet_file',type=str, help = 'file contains coordinates of outlet points')
parser_c.add_argument('-out_name',type=str, help = 'name of output solidfile (required)')
parser_c.add_argument('-dx',type=int, help = 'spatial resolution of solidfile (optional). Default is 1000')
parser_c.add_argument('-dz',type=int, help = 'lateral resolution of solidfile (optional). Default is 1000')
parser_c.add_argument('-printmask',type=int, help = 'print mask (optional). Default is 0')
parser_c.add_argument('-printbbox',type=int, help = 'print bounding box (optional). Default is 0')
#parser_c.add_argument('-z_bottom',type=int, help = 'bottom of domain (optional). Default is 0')
#parser_c.add_argument('-z_top',type=int, help = 'top of domain (optional). Default is 1000')

###required raster files
conus_pf_1k_mask = '../CONUS1_inputs/conus_1km_PFmask2.tif'
conus_pf_1k_mask = '../CONUS1_inputs/Domain_Blank_Mask.tif'

avra_path_tif = '/iplant/home/shared/avra/CONUS2.0/Inputs/domain/'
avra_path_tif = '/iplant/home/shared/avra/CONUS_1.0/SteadyState_Final/Other_Domain_Files/'

###check if file exits, if not we need to login to avra and download. This part requires icommand authorization
if not os.path.isfile(conus_pf_1k_mask):
Expand All @@ -160,8 +137,7 @@ def subset(arr,mask_arr,ds_ref, crop_to_domain, ndata=0):
if auth != 0:
print('Authentication failed...exit')
sys.exit()

os.system('iget -K '+avra_path_tif+conus_pf_1k_mask+' ../CONUS1_inputs/')
os.system('iget -vK '+avra_path_tif+tif_file+' ../CONUS1_inputs/')

###read domain raster
ds_ref = gdal.Open(conus_pf_1k_mask)
Expand Down Expand Up @@ -266,59 +242,6 @@ def subset(arr,mask_arr,ds_ref, crop_to_domain, ndata=0):
###mask array
mask_arr = (shp_raster_arr == basin_id).astype(np.int)

elif args.type == 'mask':
mask_file = args.mask_file
if not os.path.isfile(mask_file):
print (mask_file+' does not exits...please create one')
sys.exit()

file_ext = os.path.splitext(os.path.basename(mask_file))[1]
if file_ext == '.tif':
ds_mask = gdal.Open(mask_file)

#check if mask file has the same projection and extent with the domain mask file
if any([ds_ref.GetProjection() != ds_mask.GetProjection(),
np.allclose(np.array(ds_ref.GetGeoTransform()),
np.array(ds_mask.GetGeoTransform()),atol=1e-5)==False]):
print ('mask and domain do not match...exit')
sys.exit()

mask_arr = read_from_file(mask_file)

elif args.type == 'define_watershed':
dir_file = args.dir_file
if not os.path.isfile(dir_file):
print(dir_file+' does not exits...downloading from avra')
auth = os.system('iinit')
if auth != 0:
print('Authentication failed...exit')
sys.exit()

avra_path_direction = '/iplant/home/shared/avra/CONUS2.0/Inputs/Topography/Str5Ep0/'
os.system('iget -K '+avra_path_direction+dir_file+' .')

file_ext = os.path.splitext(os.path.basename(dir_file))[1]
if file_ext == '.tif':
ds_dir = gdal.Open(dir_file)

#check if direction file has the same projection and extent with the domain mask file
if any([ds_ref.GetProjection() != ds_dir.GetProjection(),
sorted(ds_ref.GetGeoTransform()) != sorted(ds_dir.GetGeoTransform())]):
print ('direction and domain do not match...exit')
sys.exit()

outlet_file = args.outlet_file
if not os.path.isfile(outlet_file):
print (outlet_file+' does not exits...please create one')
sys.exit()

dir_arr = read_from_file(dir_file)
queue = np.loadtxt(outlet_file)
queue = queue.reshape(-1,2)

#get the mask array from DelinWatershed function
mask_arr = DelinWatershed(queue, dir_arr,printflag=True)

###crop to get a tighter mask
clip_arr, new_geom, new_mask_x, bbox = subset(arr_in,mask_arr,ds_ref,crop_to_domain)

Expand Down
117 changes: 66 additions & 51 deletions Create_Subdomain/README.md
View file
Original file line number Diff line number Diff line change
@@ -1,80 +1,95 @@
# Create Subdomain
# ParFlow Mask Utilities

Create a solid file of masked domain for ParFlow
These utilities have been moved to the main Parflow repository and built when parflow is built.

### Development Team
## Building

See the list of [contributors](https://github.com/orgs/hydroframe/people) who participated in this project.
Requires C++ compiler.

## Usage
make

### Prerequisites
## Testing

To run this project on your local machine, you will need:
1. Since the code will download all the pre-processed files for the CONUS2.0 domain from Cyverse, you will need to:
* Sign up for Cyverse account
* Setting up icommands, more detail can be found [here](https://wiki.cyverse.org/wiki/display/DS/Setting+Up+iCommands)
2. Python3 installed with these specificed packages: numpy, gdal, ogr, osr, and shapely
3. *pfio* tool for reading and writing .pfb files. To download and install this tool please use this [link](https://github.com/hydroframe/tools/tree/master/pfio)
make test

### Synopsis
## mask-to-pfsol

```
python3 subset_domain.py {shapefile|mask|define_watershed} [-shp_file] [-id] [-out_name] [-dx] [-dz] [-printmask] [-printbbox] [-mask_file] [-dir_file] [-outlet_file]
```
Utitility to build 3D PFSOL domain from 2D mask file(s).

### Description
### Usage

**{shapefile|mask|define_watershed}** (Required) You need to choose one of the three subsetting methods:
1. Using a shapefile with an ID of the selected feature.
2. Using a mask file of the domain with value 1 inside the domain and value 0 elsewhere.
3. Using an output mask file from *Define_Watershed.py*. *Define_Watershed.py* defines define the watershed for a point or set of outlet points based on the flow direction file.
There are two modes of running. In the first case a single mask file
is supplied, enabling a top surface to have multiple patches with the
sides/bottom labeled as one patch. In the second mode, a mask file is
provived for each direction to enable labeling of each cell surface.
The values in the mask file are used to label the patches.

**-shp_file** (Used conjunctionally with *shapefile*) Name of the shapefile in *.shp* format. Please note that the function also requires other files in *.dbf* and *.prj* along with *.shp* file.
mask-to-pfsol
--mask <top mask input filename>
--vtk <VTK output filename>
--pfsol <PFSOL output filename>
--bottom-patch-label <bottom patch id>
--side-patch-label <side_patch id>

**-id** (Used conjunctionally with *shapefile*) ID of the selected feature within the shapefile.
Creates a PFSOL file based on 2D mask input that defines the domain
and patches on the top surface. The domain is extruded in Z to form a
3D domain. The mask input file can be many of the standard file types
supported by ParFlow, such as ParFlow binary or simple ASCI. The ASC
file format is also supported.

**-out_name** (Optional for *shapefile* and *mask*, required for *define_watershed*) Name of the output from the function. If not defined when using *shapefile* or *mask*, *out_name* will be given from selected *id* (for *shapefile*) or from *mask name* (for *mask*). When using *define_watershed*, *out_name* needed to be defined.
The mask input must be 2D with number of points in Z = 1;

**-dx** (Optional) Spatial resolution of the subset file. Default value is 1000.
The bottom and side patches are labeled with the supplied patch id's.
The top mask file is used to label patches on the top surface based on
the values in the mask file.

**-dz** (Optional) Vertical resolution of the subset file. Default value is 1000.
mask-to-pfsol
--mask-top <top mask filename>
--mask-bottom <bottom mask filename>
--mask-left <left mask filename>
--mask-right <right mask filename>
--mask-front <front mask filename>
--mask-back <back mask filename>
--vtk <VTK output filename> --pfsol <PFSOL output filename>

**-printmask** (Optional) Output an image of domain position. Default value is 0.
Each of the mask values is used to label the external boundary patches
based on the value. This enables the sides and bottom to have a more
complex patch labeling. Faces are along each axis:

**-printbbox** (Optional) Output a file which contains the bounding box of the domain. Default value is 0.
Top +Z
Bottom -Z

**-mask_file** (Used conjunctionally with *mask*) Name of the mask file. Please note that the mask file must have same extent and projection with the *input_file*.
Right +X
Left -X

**-dir_file** (Used conjunctionally with *define_watershed*) Name of the direction file. Please note that the mask file must have same extent and projection with the *input_file*.
Front -Y
Back +Y

**-outlet_file** (Used conjunctionally with *define_watershed*) Name of the outlet file. The file contains coordinates of each point (y x) by each line.
Each mask file should be the same dimensions in X and Y and have
number of points in Z = 1.

### Examples
Subsetting by mask
### ASC file format

```
python3 subset_domain.py mask -mask_file ../../mask/Coast1_mask.tif
```
The input mask is an ASC file format with the following format:

ncols 4
nrows 4
xllcorner 0.0
yllcorner 0.0
cellsize 1.0
NODATA_value 0.0
<ncols * nrows values>

Subsetting by shapefile
```
python3 subset_domain.py shapefile -shp_file ../../shp/Regions.shp -id 14
```
A 0 value is outside the domain, any other value is inside the domain.

Subsetting by a delineated watershed
```
python3 subset_domain.py define_watershed -dir_file ../../mask/Str5Ep0_direction.tif -outlet_file outlet.txt -out_name Lake2 -printmask 1
```
## pfsol-to-vtk

This utility is used to convert a PFSOL file to a VTK for easier visualization.

### Usage

pfsol-to-vtk <PFSOL input filename> <VTK output filename>


### Python workflow

A note about workflow for the python file (subset_domain.py):
1. Download all the required domain rasters from Cyverse
2. Read the mask file
3. Crop all the rasters to extents which only contain the target basin
4. Create borders (back, front, left, right, top, bottom) ascii files. **Notes:** this part adopted from Prof. Condon R codes.
5. Create a solid domain file by calling mask-to-pfsol function

Loading

0 comments on commit abcdf3c

Please sign in to comment.