2

I am attempting to create multiple watershed delineations for stations within a CSV file. I have files for flow direction (raster with a cell size of 30m), flow accumulation(raster with a cell size of 30m), and stream networks (line vector) and have snapped my points to the stream networks. Then visually checked that they do land exactly on cells with very high accumulation. Then I exported the points to a CSV file with the first column being Station ID's and the 3rd/4th column long/lat coordinates.

From there I created a loop as per this thread. but my results are coming out funny. (Screen shot below)

Each one of those black squares is not a single cell. It's 100's of columns/rows all with non-values. The extent is nowhere near as large as the delineated watershed was supposed to be so I don't think the actual delineation is within these squares.

I can't change the color of these output squares but here are an example of three of my "delineated watersheds" Screen Shot of Results

This is my python script - I run it in GrassGIS and open the results in QGIS:

import os
import sys

#set up GRASS environment variables
os.environ['GISBASE'] = "C:\\Program Files\\GRASS GIS 7.4.0"
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))

import grass.script as g
import grass.script.setup as gsetup
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
gisbase = os.environ['GISBASE']
gisdb = "D:\\nicol\\Documents\\grassdata"
location = 'newLocationNAD83'
mapset = 'PERMANENT'
gsetup.init(gisbase, gisdb, location, mapset)

g.run_command("r.in.gdal",
                input="D:\nicol\Documents\grassdata\\flow_dir_NAD83",
                output="Test_WS",
                flags="oe")

g.run_command("g.region", rast = "Test_WS", flags = "p")

#csv file reading and importation
import csv
rows = list(open("D:\\nicol\\Documents\\grassdata\\WatershedDelineation.csv"))
totalrows = len(rows) - 1

#loop through the csv(coordinates) file in r.water.outlet module
#i = start row
#j = start column
f = open("D:\\nicol\\Documents\\grassdata\\WatershedDelineation.csv")
element = list(csv.reader(f))
i = 1
j = 2
flow = "Test_WS"
basin= "WS_Station_"

while True:
    if i <= totalrows:

        g.run_command('r.wateroutlet.lessmem', 
                        input = flow , output = basin + str(element[i][j - 2]), 
                        coordinates = element[i][j+1]+","+element[i][j] )

        g.run_command('r.out.gdal', 
                        input = basin + str(element[i][j - 2]), 
                        output= "D:\\nicol\\Documents\\USA_Watersheds\\WS_Station_"+ str(element[i][j - 2])+".tif",
                        format= "GTiff",
                        type= "Float32" )
        print "Station "+str(element[i][j-2])+" Complete"
        i = i + 1

    else:
        break
  • Is it possible you have mixed up the order of the x,y coordinates in your call to r.wateroutlet ? You mentioned that column 3 is longitude, and 4 is latitude. So the coordinates= parameter should be element[i][2], element[i][3] . – Micha May 31 '18 at 11:27
  • Also, it's not clear why you are doing r.in.gdal at the beginning. Isn't the flow_dir_NAD83 already a GRASS flow direction raster? – Micha May 31 '18 at 11:29
  • The flow_dir_NAD83 is a tif file in a different location. After the first time I tried running this (and failed), I deleted those lines. I didn't have to import it again since I had "Test_WS". I quickly retyped it into the code on this form for clarification. – Nicole May 31 '18 at 15:22
  • I thought it might've been a coordinate issue. I was confused whether to use lat/long coordinates in a WGS projection or x,y coordinates in a NAD projection. Both methods output these large black squares. (If I switch the elements to element[i][2], element[i][3] then I get a single pixel that is not on this map (it's somewhere in Africa!) – Nicole May 31 '18 at 15:28
  • If this is still not working, then it's time to ask for more details. I suggest you try the following commands directly in GRASS without python. Could you post the outputs of: g.region -p, g.proj -p, r.info Test_WS and the first few lines of your CSV file WatershedDelinieation.csv ?? – Micha Jun 4 '18 at 6:18
0

The coordinates of outlets should be the same as those of the current mapset. If the current mapset is projected (not in lat/long), you need to reproject lat/long from CSV to meters or feet depending on the mapset projection using m.proj. For example, m.proj -i coordinates=long,lat and use the output coordinates in r.wateroutlet.lessmem.

Also, make sure that your flow direction map uses the same drainage encoding used by r.watershed. r.water.outlet is not compatible with ArcGIS flow directions. Please check this article.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.