mdbartos/pysheds

:earth_americas: Simple and fast watershed delineation in python.

pysheds's Issues

Catchment function doesn't check for out-of-bounds pour point

affine transformation matrix required

Hello,

Trying to follow this example example

But I fail to load the Hydrosheds data. I get a TypeError error related to the requirement of an affine transformation. The full output of the error is below. Note I am using rasterio v.0.36.0

Thanks!

 FutureWarning: The value of this property will change in version 1.0. Please see https://github.com/mapbox/rasterio/issues/86 for details.
  newinstance.read_raster(path, data_name, **kwargs)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-95-3f105dd0062f> in <module>()
----> 1 grid = Grid.from_raster(dem_data, data_name='dem')

~/dev/cust_python/lib/python3.6/site-packages/pysheds/grid.py in from_raster(cls, path, data_name, **kwargs)
    293     def from_raster(cls, path, data_name, **kwargs):
    294         newinstance = cls()
--> 295         newinstance.read_raster(path, data_name, **kwargs)
    296         return newinstance
    297

~/dev/cust_python/lib/python3.6/site-packages/pysheds/grid.py in read_raster(self, data, data_name, band, window, window_crs, metadata, **kwargs)
    282             nodata = data.dtype.type(nodata)
    283         self.add_gridded_data(data=data, data_name=data_name, affine=affine, shape=shape,
--> 284                               crs=crs, nodata=nodata, metadata=metadata)
    285
    286     @classmethod

~/dev/cust_python/lib/python3.6/site-packages/pysheds/grid.py in add_gridded_data(self, data, data_name, affine, shape, crs, nodata, mask, metadata)
    164                 pass
    165             else:
--> 166                 raise TypeError('affine transformation matrix required')
    167
    168             # initialize instance metadata

TypeError: affine transformation matrix required

TypeError: affine transformation matrix required

1. Dowloaded and unzipped this file:
   http://earlywarning.usgs.gov/hydrodata/sa_con_3s_grid/CA/n30w100_con_grid.zip

2. grid = Grid.from_raster('n30w100_con_grid/n30w100_con/n30w100_con/', data_name='dem')

3. Error: TypeError: affine transformation matrix required

Remove support for irregular grids in Grid

Network extraction doesn't remove segments contained inside other segments

Clarification: flow distance units

Hi,

Some questions about the flow_distance function.
This calculates the number of cells on the flow path from each cell in a subcatchment to the outlet of that subcatchment correct? So would multiplying a flow distance in cells by the pixel width in meters of my raster cells give a fair approximation of the flow distance in meters? Is this what specifying the weights does?

Also, in the example you give, flow_distance takes grid.catch as its input, but in the doc string it says it wants grid.dir. Which is correct?

I'm eventually planning on using this to calculate subcatchment "width" in SWMM.

Thanks!

PS: These tools are super useful thank you!

rasterio version requirements

As of v1.0a1 of rasterio, the "affine" alias for the "transform" attribute has been removed. Running Grid.from_raster gives an error DatasetReader has no attribute 'affine'. Replacing f.affine with f.transform on line 261 of grid.py seems to fix it. I don't yet know if there are other versioning issues but it would be helpful to know what version of rasterio you are using to develop the package.

Add stream order to river network extraction

Support D-infinity routing

More efficient depression fill

scikit-image provides a much more efficient way to fill depressions (nondraining flats):

http://scikit-image.org/docs/stable/auto_examples/features_detection/plot_holes_and_peaks.html#sphx-glr-auto-examples-features-detection-plot-holes-and-peaks-py

I've tested this, but need to implement it.

Implement HAND

Implement height-above nearest drainage as specified in:
https://www.sciencedirect.com/science/article/pii/S0022169411002599

error happened when install in windows 10 64bit

I came across an error when I install pysheds in windows 10 64bit using " pip install pysheds"
here is the errors:

why?
How to fix it?

Import errors

On a fresh installation of pysheds, I get import errors for geojson, affine, and pyproj.

Might trip new users. The installation of dependencies could be automated using setuptools

incorrect dinf flow accumulation in test_grid.py

Hi @mdbartos ,

While working on the test case for transport efficiency (#73), I observed a likely error in the flow accumlation with dinf.

I wanted to place a raster cell with some retention to flow near the outlet of the catchment as defined in test_grid.py. The output of the normal accumulation is attached, with d8 in blue. The reddish cell (A = 8288.15) obviously receives flow from its darkish neighbour (A=11402.86). But what about the remainder? The marked cells have accumulated values <100.

wrong variable name (likely) in _dinf_flow_distance

In the function _dinf_flow_distance, Spyder warns of an 'undefined name'. The 3rd line in this fragment should probably changed to weights_1 = weights[1].ravel(). Unfortunately, I can currently not verify the code change.

weights_0 = weights[0].ravel()
assert(isinstance(weights[1], np.ndarray))
weights_0 = weights[1].ravel() # weights_1?!
assert(weights_0.size == startnodes.size)
assert(weights_1.size == startnodes.size) # Spyder warning: 'undefined name'

error in catchment delineation, TypeError: Iterator operand 0 dtype could not be cast from dtype('float64') to dtype('int64') according to the rule 'same_kind'

Hi @mdbartos,
I use the following command :-
dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
x,y = 664976.60,5315794.1
grid.catchment(data='dir',x=x,y=y,dirmap=dirmap,out_name='catch',nodata_in=grid.nodata,nodata_out=grid.nodata)
and get the error:-

TypeError Traceback (most recent call last)
in
1 dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
2 x,y = 664976.60,5315794.1
----> 3 grid.catchment(data='dir',x=x,y=y,dirmap=dirmap,out_name='catch',nodata_in=grid.nodata,nodata_out=grid.nodata)

~/.conda/envs/bhushanPy3/lib/python3.6/site-packages/pysheds-0.2.2-py3.6.egg/pysheds/grid.py in catchment(self, x, y, data, pour_value, out_name, dirmap, nodata_in, nodata_out, xytype, routing, recursionlimit, inplace, apply_mask, ignore_metadata, **kwargs)
856 xytype=xytype, recursionlimit=recursionlimit, inplace=inplace,
857 apply_mask=apply_mask, ignore_metadata=ignore_metadata,
--> 858 properties=properties, metadata=metadata, **kwargs)
859 elif routing.lower() == 'dinf':
860 return self._dinf_catchment(x, y, fdir=fdir, pour_value=pour_value, out_name=out_name,

~/.conda/envs/bhushanPy3/lib/python3.6/site-packages/pysheds-0.2.2-py3.6.egg/pysheds/grid.py in _d8_catchment(self, x, y, fdir, pour_value, out_name, dirmap, nodata_in, nodata_out, xytype, recursionlimit, inplace, apply_mask, ignore_metadata, properties, metadata, **kwargs)
893 # get the flattened index of the pour point
894 pour_point = np.ravel_multi_index(np.array([y, x]),
--> 895 fdir.shape)
896 # reorder direction mapping to work with select_surround_ravel()
897 r_dirmap = np.array(dirmap)[[4, 5, 6, 7, 0, 1, 2, 3]].tolist()

TypeError: Iterator operand 0 dtype could not be cast from dtype('float64') to dtype('int64') according to the rule 'same_kind'

Save watershed polygon to shapefile

I found pysheds very useful for me. Thank you for sharing it on Github.
The examples are great in showing how to perform a variety of watershed processing and save the results as figures. However, showing how to save them into raster and vector files is also important. For example, after delineating a watershed, how can I save the polygon to a shapefile?

Generalize inputs for grid.clip_to

Allow the following inputs to grid.clip_to:

• Rasters
• Shapely objects

wrong variable name in _set_dirmap

Hi @mdbartos,
In _set_dirmap() you refer to an unexisting variable name direction_name in raise KeyError("{0} not found in grid instance".format(direction_name)). Would it be data, instead?

feature request: transport efficiency in flow accumulation

Hi, I propose to slightly extend the flow accumulation to consider transport efficiency, e.g. the effect of reservoirs or riparian buffers. This can be achieved by applying a correction factor between 0 and 1 in np.add.at() , e.g. np.add.at(acc, endnodes, acc[startnodes] * efficiency[startnodes]) in _d8_accumulation. Nodata would be replaced by 1 beforehand. The new argument of accumulation could either be a ndarray (like weights) or Raster/str (like flow direction). Cheers

quickstart notebook not reproducing same output

Hello, just wanted to say I think pysheds is a great project and I'm looking forward to seeing what else comes from it.

So my issue is that I am simply trying to run the quickstart notebook and when I get to the "Delineate catchment" portion I am not producing the same output. I suspect it's something to do with the pour point but I am not sure.

This is my output.

I know this is not a very informative question but I did not change anything from the quickstart notebook so I do not have much to add.

Thanks, and let me know if you want more detail.

AttributeError: 'Grid' object has no attribute 'detect_nondraining_flats'

Hi @mdbartos

I carried out "messy_dem" recipe.

depressions = grid.detect_nondraining_flats('dem')

AttributeError Traceback (most recent call last)
in
1 # Detect pits
----> 2 depressions = grid.detect_nondraining_flats('dem')

AttributeError: 'Grid' object has no attribute 'detect_nondraining_flats'

Is it possible to get array of streams in the sequence they get filled?

Streams are normally filled in order depending upon their relative hight. Many hydrological modeling number streams while creating subbasin files like BASINS etc. How can such an array be obtained in this case?!!!

Delineation of sewersheds

Great job!

Could you help on the delineation of sewersheds (subcatchments) which can be used by SWMM? Thanks!!!

Add method for filling pits and sinks in DEM

Creating flow direction and accumulation rasters using pysheds.grid methods

We are looking to create our own flow direction and accumulation rasters rather than use the coarse flow direction datasets available from HYDROSHEDS. Attached are examples of the raster I have been trying to use, but have not been successful using the flowdir function in pysheds.grids, and an example in creating these rasters would be great! I am relatively new to python and programming, so please, be gentle. Thanks! Also- are rasters saved in .tif format supported with this tool?

roi_10m.zip

Feature suggestion: burning features into DEM

Cool watershed delineation tool!

I was wondering if you're planning to add tools to manipulate DEM in the future. I think it would be neat if you could burn culvert (polyline) features or raise building elevations (polygon) into the DEM. I would be curious to see the delineation after modifications on the DEM.

dtype problem with polygonize

Using polygonize on the output of grid.catchment raises an error:

~/.conda/envs/raven/lib/python3.6/site-packages/rasterio-1.0.1-py3.6-linux-x86_64.egg/rasterio/features.py in shapes(source, mask, connectivity, transform)
    106     """
    107     transform = guard_transform(transform)
--> 108     for s, v in _shapes(source, mask, connectivity, transform):
    109         yield s, v
    110 

rasterio/_features.pyx in _shapes()

ValueError: image dtype must be one of: int16, int32, uint8, uint16, float32

grid.catchment returns an int64. Converting to int32 apparently solves the issue.

Allow raster "stitching"

accept raster / array of pour points as input to catchment

My use case is for delineating small catchments, forcing the catchments to match existing streams.
I can fill/inflate the DEM, burn stream lines into the DEM, and snap pour points to high accumulation cells. This works very well for most sites.

The exception is sites where the DEM isn't quite clean enough - there is a bump somewhere along a known stream so some catchment area is missed. Simply expanding the stream lines burned into the DEM didn't help (stream_raster = maximum_filter(stream_raster, 2)). Filtering or rounding the DEM is worth trying, or I could even burn the z-values of my known streams into the DEM. But I think the simplest way to ensure all required catchment area is captured for my use case is to:

1. generate multiple pour points along the stream
2. run the catchment function at each point
3. aggregate the results

I'll do this - but wonder if it would be worth tweaking the catchment function to also accept pour points as either a raster / array or list of points?

Windowed reading no longer works with rasterio v1.0+

Bug in check_cycles

Check cycles currently uses += n. This should be += 1 or =n

Add unit tests for raster I/O

[error] flow accumulation along nodata values differs to GRASS GIS

I test pysheds for the first time, so please be patient with me.

I have started with a flow direction grid (D8) created by the GRASS GIS tool r.watershed. The direction codes range from 1 == NE to 8 == E (actually, it also contains the negative values to indicate edge contamination, c.f. ) . The difference between the flow accumulation is zero (light blue in the attachment), as expected, except for raster cells next to nodata values (white color).

These nodata cells are necessary to exclude a channel crossing the catchment. As some streams cross the channel, I had to keep some of its raster cells, though. Any help is appreciated.

grid = Grid.from_raster(fdir, data_name='dir')
dirmap=(2,1,8,7,6,5,4,3) # coding r.watershed from N->NW
grid.dir = np.abs(grid.dir) # ignore flag edge contamination
grid.accumulation(data="dir", out_name="acc", dirmap=dirmap)

Error in messy_dem recipy: "Kernel Restarting/The kernel appears to have died. It will restart automatically."

Hi @mdbartos

in the messy_dem recipy:

grid.catchment(x, y, data='dir', out_name='catch',
               dirmap=dirmap, xytype='index')

After that kernel restart.

Support GeoTIFF?

I was wondering how I could use this with tif files. Should I convert them to ESRI GRID format first or is there any workaround for this? Thank you!

Error in delinate the catchment and grid.to_raster:"blockxsize exceeds raster width"

Hi @mdbartos,

How do you determine X and Y points?
I tried dozens different points, but i didn'n get true delinate catchment. Besides result of "grid.to_raster" process give error that "blockxsize/blockysize exceeds raster width".

Include pour point in accumulation output

Network extraction fails on forks with length one

RuntimeError: b'no arguments in initialization list'

Hi, I get this error with Windows 10, 64 bit, Anaconda 1.9.6, Spyder 3.3.2, Python 3.7.2, pysheds 0.2.3, pyproj 1.9.5.1. Please advise.

from pysheds.grid import Grid
Traceback (most recent call last):

  File "<ipython-input-16-c8708fc12346>", line 1, in <module>
    from pysheds.grid import Grid

  File "C:\ProgramData\Anaconda3\lib\site-packages\pysheds\grid.py", line 33, in <module>
    from pysheds.view import Raster

  File "C:\ProgramData\Anaconda3\lib\site-packages\pysheds\view.py", line 75, in <module>
    class BaseViewFinder():

  File "C:\ProgramData\Anaconda3\lib\site-packages\pysheds\view.py", line 77, in BaseViewFinder
    crs=pyproj.Proj('+init=epsg:4326'), y_coord_ix=0, x_coord_ix=1):

  File "C:\ProgramData\Anaconda3\lib\site-packages\pyproj\__init__.py", line 358, in __new__
    return _proj.Proj.__new__(self, projstring)

  File "_proj.pyx", line 84, in _proj.Proj.__cinit__

RuntimeError: b'no arguments in initialization list'

grid.resolve_flats is slow

Unexpected catchment output

Hi, this seems like a brilliant tool.

I am testing creating a catchment on a small DEM but not getting the expected result - only a single cell is returned: Raster([[2]]) The DEM is fine - there is definitely a contributing area higher than the pour point. My usage seems to be ok too, can you point me to something I might be missing?

Code: https://gist.github.com/smnorris/b56e8befa55cdabcd92a275bd63e79b5
DEM: https://hillcrestgeo.ca/outgoing/public/dem.tif

thanks

Allow weighted accumulation and weighted distances

Look into numba stencils

https://numba.pydata.org/numba-doc/latest/user/stencil.html

DEM Projections

I've been delineating watersheds at my job for two years now using ArcGIS but would like to switch to Pysheds. Our method involves projecting the raw DEM to an Albers Equal Area projection from geographical coordinate system first, then continuing with fill, flow direction, etc. This is "the way we do it here" but I know Pysheds doesn't involve projections before raster processing. From what I've done and read, it doesn't seem to matter much whether projection occurs before or after. Would you mind describing your reasoning for not projecting so that I can understand the process better? Does it actually matter or make a difference? I just really want to understand the process as thoroughly as I can. Thank you in advance!

[Question] Numpy derived memory error in Pysheds

Does anyone experience Memory Error when using pysheds? I am trying to follow @mdbartos example and just recreate his results to make sure things are running as-planned prior to executing the module on my own DEMs.

In this example, I can bring in the DEM and display it, but once trying to resolve_flats, that's where my Memory Error occurs. I've also recreated this error with detect_pits. It's tracing back to a function in grid.py named _inside_indices:

a = np.arange(data.size)

when data is a 6000x6000 raster, e.g. trying to develop a 36,000,000 length numpy range. The line is from grid.py, line 2010.

Anyone have any tips? My machine is a Windows 10 machine with 16gb RAM, can provide additional specs if necessary. Would code modification to include sparse arrays or looping be helpful here? If you're reading this @mdbartos , can you provide the specs of the machine you're successfully executing this code on?

Thanks!

Bug in grid.set_bbox

Recursive function d8_catchment_search causes Python to crash

Hi,

For larger DEM, the catchment method will crash Python. The cause seems to be a out-of-memory problem in the recursive function d8_catchment_search. Setting the recursion limit higher will not help and setting the recursion limit lower might avoid a hard crash, but might give a wrong (incomplete) catchment basin.

The solution would be to refactor the function from recursive to iterative. @mdbartos, I could submit an example of code when it crashes and a pull request that fixes the problem if you would guide me through the process.

Dimension error in resolve flats

The following line shouldn't be in resolve_flats:

dem_mask = np.where(dem.ravel() == nodata_in)[0]

Out of bounds error with river network extract

River network extraction yields out of bounds error when catchment is same size as view. Need to add padding to edges.

IndexError: column index (27531) out of bounds