Your Shell Continues to Use the History for the Old Mapset Grass
This page provides an overview of the steps required to pre-process and delineate a sub-basin watershed using the GRASS GIS python shell. The test domain used in this example is the Little Pic River Near Coldwell (WSC station 02BA003).
This workflow follows the same steps as the workflow using the GRASS GUI, but uses the python shell within GRASS to complete the processing. The scripting file can be found here.
Opening the Python shell within GRASS GIS.
Selecting the python logo in the GRASS Layer Manager opens the GRASS GIS Python Editor
Python Script
Once the GRASS python shell is open, copy and paste the script below into the editor and selectrun. The progress of the script will be displayed in theconsoletab, and any errors that occur will be listed here.
Setting paths and downloading the required GRASS modules
The following code loads the necessary libraries required by the code and sets the paths to the correct input and output data sources.
Setting arguments and importing DEM and WSC data
The code below outlines what arguments are required for the script. The threshold refers to theminimum flow accumulation used to generate the stream network. The user must also input the correct WSC basin ID and the coordinates of the outlet.
Discretion must be used when deciding what minimum flow accumulation value to use, as it will impact the number of stream and sub-basins produced by the discretization.
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE #Required arguments WSC_basin = 'EC_02BA003_1' #Water Survey Canada (WSC) Basin ID coords = '-86.60734847236674,48.84927450473443' #Coordinates of the outlet threshold = '8500' #Minimum flow accumulation for streams # Import Raster DEM_in = 'DEM_in' gcore.run_command('r.in.gdal', input = input_DEM, output = DEM_in, overwrite = True) #Import geodatabase Basin_in ='Basin_in' gcore.run_command('v.import',layer = WSC_basin, input = input_WSC, output = Basin_in, overwrite = True) Buffering basin and producing the watershed outline
The following script generates a buffered outline of the WSC basin as well as a watershed outline using theflow directionand outlet coordinates.
At two points in this script a MASK is set, which willblock out certain areas of a raster map from analysis and/or display, by "hiding" them from sight of other GRASS modules. Data falling within the boundaries of the MASK can be modified and operated upon by other GRASS raster modules; data falling outside the MASK is treated as if it were NULL.
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE #buffer basin basin_buffered = 'basin_buffered' gcore.run_command('v.buffer', input = Basin_in, distance = '0.05', minordistance = '0.05', output = basin_buffered, overwrite = True) #update mask #add if mask true, then remove gcore.run_command('r.mask', vect=basin_buffered, overwrite = True) #set computational region #Note: the next step will fail if the computational region is not set correctly gcore.run_command('g.region', raster = DEM_in, vector = basin_buffered) #generate flow accumulation and direction flow_dir = 'flow_dir' flow_acc = 'flow_acc' gcore.run_command('r.watershed', elevation = DEM_in, accumulation = flow_acc, drainage = flow_dir, flags = 's', overwrite = True) #generate watershed outline watershed_outline ='watershed_outline' gcore.run_command('r.water.outlet', input = flow_dir, output = watershed_outline, coordinates = coords, overwrite = True) #update mask gcore.run_command('r.mask', rast=watershed_outline, overwrite = True) Generating stream network and topology (Rank, Next) and generating sub-basins
The following script generates the stream network and topology and sub-basins.The flow accumulation raster may be used to estimate the m inimum flow accumulation for stream. In this example, a minimum flow accumulation value of9016.44 provides a reasonable discretization of thestream network. r.stream.segment is used to generate the stream toplogy information (RANK, NEXT).The module r.stream.basins is used to delineate basins and sub-basins with different input data. In this example, an input stream raster and flow direction will be used to delineate the watershed into sub-basins. Refer to the r.stream.basins documentationhere for more examples and use cases.
Bothr.stream.segmentandr.stream.basinsare both add-on modules.g.extensionis run at the beginning of the script to install the modules.
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE #Generate stream network stream_raster = 'stream_raster' gcore.run_command('r.stream.extract', elevation = DEM_in, accumulation = flow_acc, threshold = threshold, stream_raster = stream_raster, direction = flow_dir, overwrite = True) #generate stream topology stream_segments = 'stream_segments' stream_sector = 'stream_sector' gcore.run_command('r.stream.segment', stream_rast = stream_raster, elevation = DEM_in, direction = flow_dir, segments = stream_segments, sectors = stream_sector, overwrite = True) #Generate subbasins watershed_subbasins ='watershed_subbasins' gcore.run_command('r.stream.basins', direction = flow_dir, stream_rast = stream_raster, basins = watershed_subbasins, overwrite = True) Cleaning the data and converting to a vector dataset
Sometimes r.stream.basins module will create sub-basins that are not one continuous feature, such that a small piece may be 'outside' of the larger sub-basin. The following script uses the r .neighbors and gdal.sieve toolsprior to converting watershed_subbasin to a vector . In some situations, the r.neighbors tools is sub-optimal and will not adequately 'clean' the raster. Because of this, the gdal.sieve tool is used as well. The gdal.sieve script removes raster polygons smaller than a provided threshold size (in pixels) and replaces them with the pixel value of the largest neighbour polygon. In this example, a threshold of 10 is used, but this value can be adjusted if needed.
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE #Clean data watershed_subbasins_cleaned = 'watershed_subbasins_cleaned' gcore.run_command('r.neighbors', input = watershed_subbasins, output = watershed_subbasins_cleaned, method = 'mode', overwrite = True) ########################################################### # sometimes r.neighbors does not remove all spilt polygons, the tool below (gdal.sieve) completes a similar process. # the threshold parameter can be adjusted based on the DEM. #Export watershed_subbasins gcore.run_command('r.out.gdal', input = watershed_subbasins_cleaned, output = proj_path +'watershed_subbasins_cleaned.tif', format = 'GTiff', overwrite = True) #import raster as gdal object Image = gdal.Open(proj_path + 'watershed_subbasins_cleaned.tif', 1) # open image in read-write mode Band = Image.GetRasterBand(1) gdal.SieveFilter(srcBand = Band, maskBand = Band, dstBand = Band, threshold = 10) #Export gdal.GetDriverByName('Gtiff').CreateCopy(proj_path + 'watershed_subbasins_sieved.tif', Image, strict = 0) #import the sieved file watershed_subbasins_sieved = 'watershed_subbasins_sieved' gcore.run_command('r.in.gdal', input = proj_path + 'watershed_subbasins_sieved.tif', output = watershed_subbasins_sieved, overwrite = True) ################################################################ #Convert to vector subbasins_vector = 'subbasins_vector' gcore.run_command('r.to.vect', input = watershed_subbasins_sieved, output = subbasins_vector, type ='area', flags = 'v', overwrite = True) Joiningstream_segmentsandsubbasin_vector and calculating Area, Latitude and Longitude
The script below adds the required columns to the database (area, lat, lon), and populates the attribute table using the v.to.dbmodule.
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE gcore.run_command('v.db.join', map = subbasins_vector, column='cat', other_table =stream_segments, other_column = 's_order', overwrite = True) #add new columns for area, lat and lon gcore.run_command('v.db.addcolumn', map = subbasins_vector, columns = "area double precision") gcore.run_command('v.db.addcolumn', map = subbasins_vector, columns = "lat double precision") gcore.run_command('v.db.addcolumn', map = subbasins_vector, columns = "lon double precision") #generate values for new columns gcore.run_command('v.to.db', map = subbasins_vector, option = 'area', columns = 'area') gcore.run_command('v.to.db', map = subbasins_vector, option = 'coor', columns = 'lon,lat') Exporting the attribute table as a.csv
If no landcover information is required in the drainage database, then the final step is to export the attribute table as a .csv file
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE #Export attribute table gcore.run_command('db.out.ogr', input = subbasins_vector, output = output_db, overwrite = True) (Optional) Processing CLASS landcover data
The landcover raster must be addedbefore running the grass script below. Refer to the instructions here.
Importing landcover raster from thefirstlocation
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE #Define output file output_LC_db = r'C:\Users\DisherB\Documents\Watershed_Delin\Python-Dev\drainage_LC_out.csv' #import landcover raster from location #note: the due to an error with raster reprojection in GRASS, the landcover raster must first be added to a new location #and imported into this session. Refer here for instructions: #location = location containing LC data #mapset = mapset name (usually PERMANENT) #input = name of landcover data location = 'Landcover' mapset = 'PERMANENT' input = 'CAN_NALCMS_2015_v2_land_cover_30m' gcore.run_command('r.proj', location = location, mapset = mapset, input = input, overwrite = True) Reclassifying and converting the raster to vector
MESH requires that the number of classes for the CLASS LSS be reduced. Typically, the number of classes is reduced to seven (following the table here), but may also be customized by the user. The fileCLASS_reclass_rules.txtcan be found within GitHub repository.
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE #reclassify raster data #Note: the file used in this example is found within the github file_path = r'C:\Users\DisherB\Documents\Watershed_Delin\Watershed_delineation\Code\Network Topology\Input\Bow\CLASS_reclass_rules.txt' gcore.run_command('r.reclass', input = input, output = 'Landcover_reclass', rules = file_path, overwrite = True) #convert to vector gcore.run_command('r.to.vect', input = 'Landcover_reclass', output = 'Landcover_reclass_vec', type = 'area', overwrite = True) Overlaying and calculating area
Next, the reclassified landcover classification must be overlaid onto basin file created in the previous page, and the area must be determined for each for each parcel in the attribute table.
# THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE # THIS SCRIPT CONTINUES FROM ABOVE #Overlay sub-basin delineation and landcover data gcore.run_command('v.overlay', ainput = 'subbasins_vector', binput = 'Landcover_reclass_vec', operator = 'and', output = 'LC_overlay', overwrite = True) #add new columns for area gcore.run_command('v.db.addcolumn', map = 'LC_overlay', columns = "area double precision") #generate values for new columns gcore.run_command('v.to.db', map = 'LC_overlay', option = 'area', columns = 'area') Exporting the attribute table
The final step in GRASS is to export the attribute table as a .csv file. This file will be processed in python to complete the final MESH drainage file.
#Export attribute table gcore.run_command('db.out.ogr', input = 'LC_overlay', output = output_LC_db, overwrite = True) Source: https://wiki.usask.ca/display/MESH/Using+the+GRASS+python+shell+to+pre-process+the+DEM
Belum ada Komentar untuk "Your Shell Continues to Use the History for the Old Mapset Grass"
Posting Komentar