Mapping subcatchments to nodes

landuse2inp.py

# -*- coding: utf-8 -*-
"""
@author: Sebastian Schoeller

Source:
https://github.com/Deltares/pyflwdir/issues/15

1.    Create regions from the landuse map, using e.g. scipy.ndimage.label. This returns a map with unique IDs for connected regions with the same landuse.
2.    Find the most downstream "outflow pixel" of each region using pyflwdir.FlwdirRaster.inflow_idxs Note: should still be added to the API Reference in the docs. This returns a list with linear indices of outflow pixels.
3.    Derive subbasins using pyflwdir.FlwdirRaster.basins based on the indices of the outflow pixels from step 2, supplied using the idxs_out option. This returns a map with unique IDs for each basin.
4.    Find the relations between the basins using pyflwdir.FlwdirRaster.streams, again based on the indices of the outflow pixels from step 2, supplied using the idxs_out option. This returns geofeatures with river segments between outflow pixels and per segment the index of the next downstream segment (and thus basin).

https://numpy.org/devdocs/reference/generated/numpy.ma.masked_where.html
https://numpy.org/doc/stable/reference/generated/numpy.nanmax.html
https://xarray.pydata.org/en/stable/generated/xarray.open_rasterio.html
https://automating-gis-processes.github.io/CSC18/lessons/L6/clipping-raster.html
"""

import io
import geopandas as gpd
import numpy as np
import rasterio
from rasterio.plot import show
from rasterio import features
from geocube.api.core import make_geocube
from pandas import read_csv
import pyflwdir
import rioxarray
import xarray
from xrspatial.zonal import regions
from affine import Affine
import matplotlib.pyplot as plt

lu = """id,imperv,dst_imp,n_imper,dst_per,n_perv,percz_i,rain_ga,conduct,initdef,suction,Snowpacks,Tag
10,100,0.87,0.012,0.87,0.012,0,STA01,4.21,0.217,88.9,sp_imp_nonplowable,Dach
12,100,0.87,0.012,0.87,0.012,0,STA01,4.21,0.217,88.9,sp_imp_nonplowable,Kiesdach
14,100,0.87,0.012,0.87,0.012,0,STA01,4.21,0.217,88.9,sp_imp_nonplowable,Gründach
26,100,0.87,0.012,0.87,0.012,0,STA01,4.21,0.217,88.9,sp_imp_nonplowable,Container
30,100,0.42,0.011,0.42,0.011,0,STA01,4.21,0.217,88.9,sp_imp_plowable,Asphalt
32,100,0.39,0.02,0.39,0.02,0,STA01,4.21,0.217,88.9,sp_imp_plowable,Pflaster
34,100,4.22,0.238,4.22,0.238,0,STA01,4.21,0.217,88.9,sp_imp_nonplowable,Senke
60,0,4.22,0.238,4.22,0.238,0,STA01,4.21,0.217,88.9,sp_per,Grünfläche
62,0,4.22,0.238,4.22,0.238,0,STA01,4.21,0.217,88.9,sp_per,Wald
64,100,4.22,0.238,4.22,0.238,0,STA01,4.21,0.217,88.9,sp_imp_nonplowable,Wasser
"""

# Print summary of raster information
def print_raster(raster):
    print(
        f"shape: {raster.rio.shape}\n"
        f"resolution: {raster.rio.resolution()}\n"
        f"bounds: {raster.rio.bounds()}\n"
        f"sum: {raster.sum().item()}\n"
        f"CRS: {raster.rio.crs}\n"
    )

# Reading green-ampt mapping table to dataframe
landuse = read_csv(io.StringIO(lu), index_col=0)

# Reading landuse to dataframe
gdf_lu = gpd.read_file(r'dataset.gpkg', layer='landuse')
gdf_lu.plot(column='type')
bounds = gdf_lu.total_bounds

# Reading junctions to dataframe
gdf_ju = gpd.read_file(r'dataset.gpkg', layer='junctions')
gdf_ju.plot(column='elevation')

# Reading dem from file
with rioxarray.open_rasterio(r'dataset.gpkg', layer='dem').sel(band=1) as src1:
    # elevtn = src.read(1)
    nodata = src1.rio.nodata
    transform = src1.rio.transform
    # extent = np.array(src1.rio.bounds)[[0, 2, 1, 3]]
    crs = src1.rio.crs
    latlon = crs.to_epsg() == 4326
    # prof = src.profile
    
with rasterio.open(r'dataset.gpkg', layer='dem') as src:
    elevtn = src.read(1)
    nodata = src.nodata
    transform = src.transform
    extent = np.array(src.bounds)[[0, 2, 1, 3]]
    crs = src.crs
    latlon = crs.to_epsg() == 4326
    prof = src.profile
show(elevtn)

# 1.Create regions from the landuse map
# Convert dataframe to raster
ds_lu = make_geocube(
    vector_data=gdf_lu,
    output_crs=gdf_lu.crs.to_epsg(),
    resolution=(-1, 1),
)

# 1a.Project DTM to landuse dataset sizes 
ds_dtm = src1.rio.reproject_match(ds_lu).to_dataset(name='dtm')

# 1b.Assign coordinates
ds_dtm = ds_dtm.assign_coords({
    "x": ds_lu.x,
    "y": ds_lu.y,
})

# # 1c.Compare visual
# fig, axes = plt.subplots(ncols=2, figsize=(12,4))
# ds_lu.plot(ax=axes[0])
# ds_dtm.plot(ax=axes[1])
# plt.draw()

# 1d.Create regions
ds_lu_rg = regions(raster=ds_lu.type,neighborhood=8).to_dataset()

# 1e.Unify datasets to one for further use
# ds=ds_lu.assign(elevtn=ds_dtm.band)
# ds=ds.assign(region=ds_lu_rg.regions)

# 1f.Print summary info
np.nanmax(ds_lu_rg.to_array)
np.ma.masked_where(ds_lu_rg.to_array == 614, ds_lu_rg.to_array)

# 1g.Testing raster coordinates
# Error as per https://github.com/corteva/rioxarray/issues/209 needs fixing
ds_dtm.rio.to_raster("test.tif")
ds_lu.rio.to_raster("test1.tif")
ds_lu_rg.rio.to_raster("test2.tif")

# 1f.Create pyflw array
flw = pyflwdir.from_dem(
    data=ds_dtm,
    nodata=ds_dtm.rio._nodata,
    transform=ds_dtm.rio.transform,
    latlon=ds_dtm.rio.crs.to_epsg(),
    outlets="min",
)

# Testing
# Assign elevation to dataset "ds"
ds=ds_lu.assign(elevtn=src)
# Save to geotif
#ds.rio.to_raster("test.tif")
# Create regions
lu_rg_ar = regions(raster=ds_lu.type,neighborhood=8).to_numpy()
# Return number of regions
lu_rg_no = np.nanmax(lu_rg_ar)
# Create mask
mask = np.ma.masked_where(lu_rg_ar == 614, lu_rg_ar)