Fast methods to work with hydro- and topography data in pure Python.

Dear @DirkEilander ,

This is Jiangchao Qiu, a Ph.D. student from Sun Yat-sen University, China. I am major in simulating storm surge using physical model.

I know you by reading your paper : The effect of surge on riverine flood hazard and impact in deltas globally. Great works!

Recently, my collaborator provides me some global modeling result using CaMa-Flood in 15 arcmin resolution. I am trying to analyze these result and make some visualization, especially the river network and river mouths along the coastlines.

Firstly, I used the function of cmf_maps_io.py (in your compound_hotpots repositories) to transform the nextxy.bin file to nextxy.tif file glb_15min.zip

Secondly, I want to use the pydlwdir package to visualize the river network, but some errors occurred when I parse the nextxy.tif file to
the actionable common format.

the following is the screensnaps

I try to use the function of pyflwdir.FlwdirRaster and pyflwdir.from_array, but both failed, I don't know how to fix it, could you please help me check and give some suggestion about how to solve this problem.

By the way, the environment of my package seems no problem, since all of the example in the notebook folder can be carried out smothly.

Many thanks to you and looking forward to your reply.

Best regards, Jiangchao

Hi @DirkEilander ,

I want to identify all of the river mouth and make a plot using the 3 arcmin river network dataset nextxy.bin.

But it is odd that the figure I plotted here seems reversed in the latitude. Could you please have a look? Many thanks.

Best regards, Jiangchao

my scripts:

import pyflwdir import numpy as np

bbox=[-180, -90, 180, 90] data, transform = pyflwdir.read_nextxy( fn='glb_03min/nextxy.bin', nrow=3600, ncol=7200, bbox=bbox)

import xarray as xr xr.DataArray(data[0]).plot()

x = data[0] lon = np.linspace(-180,180,7200) lat = np.linspace(-90,90,3600) lon[np.where(x==-9)[1]],lat[np.where(x==-9)[0]]

import geopandas as gpd from shapely.geometry import Point gdf = gpd.GeoDataFrame([Point(x, y) for x,y in zip(lon[np.where(x==-9)[1]], lat[np.where(x==-9)[0]])]) gdf.columns = ['geometry']

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) fig, ax = plt.subplots() gdf.plot(ax=ax,markersize=0.1, color='red') world.plot(ax=ax);

My apologies for opening this issue as I don't know how else to do this, but I was wondering if the pyfldir toolkit can be used to calculate the distance to a river network? I've been desperately looking for an open-source tool to do this. So a method that uses a flow direction raster to calculate the flow distance to the nearest stream.

Feature request: This may be in the code but I haven't see it in the docs anywhere. The flow path trace function returns a list of pixel indices in the linear index format. A helper function to convert to the coordinate system of the raster would be helpful.

Dear all,

thanks for this package. I am running conda 4.11 on Win10 and have installed pyflwdir as instructed. When following the user guide instructions for plotting I stumble upon below error.

Python 3.8.12 | packaged by conda-forge | (default, Oct 12 2021, 21:22:46) [MSC v.1916 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> from utils import quickplot, colors, cm
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'utils'

Any hints on how to alleviate the situation are highly welcome. Same counts for install under Ubuntu WSL.

Kind regards

Sebastian

I think I have the error narrowed down to how the JIT compiler views either the rasterio object or the parameters sent to pyflwdir.from_dem() function. I have a simple example based on the Quickstart from the documentation. The error occurs when i try to create a flow direction raster from elevation read through rasterio.
https://colab.research.google.com/drive/1engswQ4uaKB0UMMM4zt8YzC8wEAnK1o1?usp=sharing

• FlwdirRaster.streams method includes a (zero-length) line at pits and adds option for trace direction in combination with segment end cells.
• Moved accuflux method from FlwdirRaster to parent Flwdir class.

Subgrid river lengths can locally be very short which slows down the model. Potentially we can redistribute the subgrid river lengths within a single river branch to speed up calculations without making a large computational error.

method in flwdir.py

• inputs: river length, min threshold

numba jitted method in streams.py

• input: idxs_ds; idxs_seq (sorted from down- to upstream); confluences (based on number of upstream cells), river lengths, min threshold

suggested algorithm

from up- to downstream:
   while not at next confluence or less than n cells:
      get lengths & indices
   average lengths

assert that the total river length is the same

np.random.seed(11)
dem = np.random.rand(15, 10)
flwdir = pyflwdir.from_dem(dem, outlets='min')
dem1 = flwdir.dem_adjust(dem)
np.all((dem1 - flwdir.downstream(dem1))>=0)
>> False

a method that maps subgrid cells to river cell (at model resolution) outflow points based on D8, sorts the HAND values (in future also manning?) and returns the contributing area for fixed HAND intervals

Hi,

I would like to create polygons of the contributing area for > 10,000 points. I have all information that I need for this, including a 30m DEM, a river network (LineString geodataframe file) and derived flow direction raster. All compiled using pyflwdir.

I'm currently using the pyflwdir tool, but this takes around 5 minutes per point.. Has anyone any clues on how to approach this easier and quicker? Also, to insert the streams from a LineString geodataframe (instead of this flow.stream_order() > 4 command). This is my current script:

with rasterio.open(flwdir_fn, "r") as src: flwdir = src.read(1) crs = src.crs flw = pyflwdir.from_array( flwdir, ftype="d8", transform=src.transform, latlon=crs.is_geographic, cache=True, )

for point in geo_df["geometry"].head(): x = point.x y = point.y print(x, y) subbasins = flw.basins(xy=(x, y), streams=flw.stream_order() >= 4)

Many many thanks!

Dear all,

I am trying to resemble a functionality similar to GisToSWMM using pyflwdir. In my understanding of this algorithm a subcatchment hierarchy is determined considering landuse. Every (sub-)subcatchment contains only one landuse and routes either to another (sub-)subcatchment or an outlet. Outlets are given in the form of enlarged nodes, thus rasterized sink information derived from outlet nodes and DEM sinks.

Reading through pyflwdir documentation I figure that the Pfafstetter method allows for deriving hierarchical subcatchment information. If understood correctly it permits to figure out a form of subcatchment routing, thus which subcatchment drains to another subcatchment or outlet.

I am seeking ways to mask this routing method using landuse, so that every subcatchment contains only one landuse. Any ideas if and/or how to achieve this goal are highly appreciated.

Kind regards

Sebastian

• streams array can be simplified, no need to save cell indices
• lots of duplicate code between upscale_check and upscale_error methods
• remove unused code
• improve comments. consistent use of pixel vs cell in terminology.