use pyflwdir package to parse nextxy data of CaMa-Flood!

[Jiangchao3]

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

[Jiangchao3]

And I also email to Prof. Dai Yamazaki, he replied that the traditionary GMT and Fortran code for river network visualization don't work now, so I am trying to learn the pyflwdir toolbox recently.

[DirkEilander]

Hi @Jiangchao3,

Thanks for using pyflwdir!

I had a look at your data and noticed that pits have a value of -10 instead of -9 in the nextxy data. Perhaps this has changed in the recent CaMa-Flood updates? Pinging Prof Yamazaki (@bigasmountain) here to confirm this has indeed changed in which case I'll update the code.

For now you can just change the pit values in the nextxy data to use pyflwdir.

import pyflwdir
import rasterio
import numpy as np

# read data
with rasterio.open('nextxy.tif') as src:
    data = src.read()
    transform = src.transform
    latlon = src.crs.is_geographic

# change pit values and parse with pyflwdir
data = np.where(data==-10, -9, data)
flw = pyflwdir.from_array(data, ftype='nextxy')

To make nice visualization you will want to use either the vectorize or streams method. Both methods can take an additional xs and ys argument with the outlet x and y coordinates which will make your visualizations look much better.

Hope this helps!
Dirk

PS. For a next issue anywhere on github it is always useful if you paste the actual code instead of a screenshot to quickly be able to reproduce it. Otherwise your issue was really clear and had the data for me to debug it!

[bigasmountain]

Yes -10 is inland termination, while -9 is for usual river mouth flowing to ocean,

[Jiangchao3]

Hi @DirkEilander and @bigasmountain, thanks very much for your kindly reply.
I will follow your suggestions to paste the code instead of screenshot on any of my next github issue.

And yes I also noticed that the latest CaMa-Flood manual describes that -10 represents the inland termination while -9 represents the river mouth.

[Jiangchao3]

Hi @Jiangchao3,

Thanks for using pyflwdir!

I had a look at your data and noticed that pits have a value of -10 instead of -9 in the nextxy data. Perhaps this has changed in the recent CaMa-Flood updates? Pinging Prof Yamazaki (@bigasmountain) here to confirm this has indeed changed in which case I'll update the code.

For now you can just change the pit values in the nextxy data to use pyflwdir.

import pyflwdir
import rasterio
import numpy as np

# read data
with rasterio.open('nextxy.tif') as src:
    data = src.read()
    transform = src.transform
    latlon = src.crs.is_geographic

# change pit values and parse with pyflwdir
data = np.where(data==-10, -9, data)
flw = pyflwdir.from_array(data, ftype='nextxy')

To make nice visualization you will want to use either the vectorize or streams method. Both methods can take an additional xs and ys argument with the outlet x and y coordinates which will make your visualizations look much better.

Hope this helps! Dirk

PS. For a next issue anywhere on github it is always useful if you paste the actual code instead of a screenshot to quickly be able to reproduce it. Otherwise your issue was really clear and had the data for me to debug it!

Hi @DirkEilander , your recomended script worked after replace the pits of -10 with -9. I will spend some time to further learn the whole package of pyflwdir. And looking forward to your update of the code to be able to parse the new river network dataset of cama-flood. Many thanks for your kindly help.

[DirkEilander]

This is fixed in 16b1ac0 (still unrealeased - for now you get the update by updating pyflwdir directly from git)
In addition I've added a convenience method to read the binary CaMa-Flood nextxy data.
The following example should provide a nice plot of your CaMa-Flood rivers.

# %%
import pyflwdir

# %% read data
bbox=[-180, -90, 180, 90]
data, transform = pyflwdir.read_nextxy(
    fn='nextxy.bin', nrow=720, ncol=1440, bbox=bbox
)
flw = pyflwdir.from_array(
    data,
    ftype='nextxy',
    latlon=True,
    transform = transform

)
# %% vectorize streams
import geopandas as gpd
feats = flw.streams(strord=flw.stream_order(), min_sto=2)
gdf = gpd.GeoDataFrame.from_features(feats, crs=4326)

# %% drop lines crossing anti-meridian
import numpy as np
def check_crossing(line):
    return np.any(np.abs(np.diff([xy[0] for xy in line.coords[:]]))>180)
drop_lines = np.vectorize(check_crossing)(gdf.geometry)
gdf1 = gdf[~drop_lines]

#%% plot values by stream order
gdf1.sort_values('strord').plot(column='strord')

[Jiangchao3]

So great for providing such a detailed script. I will check it after getting back to office tomorrow! Thanks again!

[Jiangchao3]

Hi @DirkEilander, the above script you provided works well, many thanks to your sincere help!