Comments (11)
import heapq
import numpy as np
def priority_flood(dem):
open = []
closed = np.zeros_like(dem, dtype='bool', subok=False)
y, x = dem.shape
open.extend(zip(dem[0], range(x), repeat(0)))
open.extend(zip(dem[y-1], range(x), repeat(y-1)))
open.extend(zip(dem[:, 0], range(y), repeat(0)))
open.extend(zip(dem[:, x-1], range(y), repeat(x-1)))
heapq.heapify(open)
closed[0] = True
closed[y-1] = True
closed[:, 0] = True
closed[:, x-1] = True
pq_pop = heapq.heappop
pq_push = heapq.heappush
row_offsets = np.array([-1, -1, 0, 1, 1, 1, 0, -1])
col_offsets = np.array([0, 1, 1, 1, 0, -1, -1, -1])
while open:
c, y1, x1 = pq_pop(open)
cols = x + col_offsets
rows = y + row_offsets
for r, c in zip(rows, cols):
if 0 <= r < y and 0 <= c < x and not dem.mask[r, c]:
if closed[r, c]:
continue
dem[r, c] = max(dem[r, c], dem[y1, x1])
closed[r, c] = True
pq_push(open, (dem[r, c], r, c))
return dem
It's a direct implementation of Algorithm 1 in the linked paper.
from pysheds.
Not sure. Did you try setting the epsilon value and number of iterations as suggested in #238?
Just making sure that you are doing fill_pits > fill_depressions > resolve_flats in that order.
I could take a shot at implementing the methods as they are programmed in richdem. Which methods specifically?
- https://github.com/r-barnes/richdem/blob/master/programs/tiled_flat_resolution/main.cpp
- https://github.com/r-barnes/richdem/blob/master/programs/parallel_priority_flood/Zhou2016pf.hpp
- https://github.com/r-barnes/richdem/blob/master/programs/parallel_priority_flood/Barnes2014pf.hpp
from pysheds.
Another potential issue is that in this case a single depression/flat is touching multiple edges (3 in this case).
from pysheds.
I should mention, that this issue arises on specific DEMs. I have observed that pysheds seems to struggle with DEMs that don't have dramatic changes in elevation. These set of DEMs that I'm working with are coastal areas.
The plots above are a 1000x1000 window of a much larger DEM (8000x8000 roughly). The smaller window was for plotting purposes only. The methods were run on the full DEM.
I am doing the operations in the correct order (except I'm not filling pits, I thought those were also handled by fill_depressions). So the order is fill_depressions > resolve_flats. However, in this specific scenario I was trying to isolate things, so the fill depressions plot is only fill_depressions and resolve flats is only resolve_flats. Though, running both fill_depressions and resolve_flats in sequence doesn't change the nature of the results. Each method in both packages were run with the default arguments. I didn't change the value of epsilon. In richdem, I am using the equivalent methods of FillDepressions > ResolveFlats.
from pysheds.
See if the methods suggested in #238 solve the issue. The fact that it is struggling with flat areas makes me think that lowering the epsilon value could help.
from pysheds.
@mdbartos I tried adjusting the epsilon value when resolving flats, however, it doesn't seem to fix the issue.
I then tried reimplementing fill_depressions
using the Priority-Flood algorithm and I now get results very similar to richdem. I think the issue is using sk-image to fill depressions.
I can put together a PR later today or tomorrow.
from pysheds.
Hello @groutr ! I had the same problem as described in this issue. Unfortunately, the tips suggested by @mdbartos didn't bring any improvement. I was wondering if you could kindly share your implementation of fill_depression
mentioned in the previous comment 🥺
from pysheds.
Hi @groutr, this is excellent. Would you be interested in creating a PR?
from pysheds.
@mdbartos The code above was simply for testing agaist pysheds. For a PR, I wanted to implement the better, but slightly more complex algorithms in the paper, but unfortunately haven't had time yet.
Ultimately, I think Alg 2 (better priority-flood) and Alg 4 might be useful for pysheds.
from pysheds.
I'm not familiar with optimizing numba code. My best attempt processes a DEM with shape of (12150, 9622) in about 40min. I suspect the heapq-based priority queue is really slowing things down.
Just for kicks, I wrote a cython version and it runs in about 10min on the same DEM.
from pysheds.
I recently saw this tool for profiling numba code: https://pythonspeed.com/articles/numba-profiling/
Could be useful for diagnosing the slowdown.
from pysheds.
Related Issues (20)
- How to create Flow Distance raster same shape as DEM raster?
- Re-evaluate Numba performance HOT 2
- Allow Numba to perform polymorphic dispatching HOT 5
- Can't 'imshow' DEM for some unkown reason HOT 1
- extract_profiles function returning wrong connections when neighbouring cells drain to outlet
- strang plots HOT 1
- Allow user to disable `parallel=True` with numba HOT 2
- Why stream network `LineStrings` do not pass through the centroid of each grid cell? HOT 1
- accumulation issue HOT 3
- Using catchment() with discharge point outside of extents crashes program
- sGrid and pGrid have major differences
- Pysheds Cupy / Cuspatial support HOT 1
- Accumulation disconnections HOT 4
- setup fail by github repo of latest release version
- Pixels that should be accumulation watercourses are shown as nodata HOT 2
- Wrong bbox when using `ViewFinder` HOT 3
- issue in pygrid with np.int and np.warnings
- D8 and Dinf flow directions look incorrect HOT 4
- Indexing Error in _d8_distance_from_ridge
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from pysheds.