• 基于Pyflwdir实现流域的提取(参照官网例子)


    本文参照官网例子实现流域的提取,官方GitHub地址如下pyflwdir

    该工具包目前仅支持D8和LDD两种算法,在效率上具有较好的应用性,我用省级的DEM(30米)数据作为测试,输出效率可以满足一般作业需要。

    环境environment.yml如下:

    1. name: pyflwdir
    2. channels:
    3. - conda-forge
    4. # note that these are the developer dependencies,
    5. # if you only want the user dependencies, see
    6. # install_requires in setup.py
    7. dependencies:
    8. # required
    9. - python>=3.9
    10. - pyflwdir
    11. # optional for notebooks
    12. - cartopy>=0.20
    13. - descartes
    14. - geopandas>=0.8
    15. - jupyter
    16. - matplotlib
    17. - rasterio

    用到的utils.py的代码如下:

    1. import matplotlib
    2. import matplotlib.pyplot as plt
    3. from matplotlib import cm, colors
    4. import cartopy.crs as ccrs
    5. import descartes
    6. import numpy as np
    7. import os
    8. import rasterio
    9. from rasterio import features
    10. import geopandas as gpd
    11. np.random.seed(seed=101)
    12. matplotlib.rcParams["savefig.bbox"] = "tight"
    13. matplotlib.rcParams["savefig.dpi"] = 256
    14. plt.style.use("seaborn-v0_8-whitegrid")
    15. # read example elevation data and derive background hillslope
    16. fn = os.path.join(os.path.dirname(__file__), r"D:\CLIPGY.tif")
    17. with rasterio.open(fn, "r") as src:
    18. elevtn = src.read(1)
    19. extent = np.array(src.bounds)[[0, 2, 1, 3]]
    20. crs = src.crs
    21. ls = matplotlib.colors.LightSource(azdeg=115, altdeg=45)
    22. hs = ls.hillshade(np.ma.masked_equal(elevtn, -9999), vert_exag=1e3)
    23. # convenience method for plotting
    24. def quickplot(
    25. gdfs=[], raster=None, hillshade=True, extent=extent, hs=hs, title="", filename=""
    26. ):
    27. fig = plt.figure(figsize=(8, 15))
    28. ax = fig.add_subplot(projection=ccrs.PlateCarree())
    29. # plot hillshade background
    30. if hillshade:
    31. ax.imshow(
    32. hs,
    33. origin="upper",
    34. extent=extent,
    35. cmap="Greys",
    36. alpha=0.3,
    37. zorder=0,
    38. )
    39. # plot geopandas GeoDataFrame
    40. for gdf, kwargs in gdfs:
    41. gdf.plot(ax=ax, **kwargs)
    42. if raster is not None:
    43. data, nodata, kwargs = raster
    44. ax.imshow(
    45. np.ma.masked_equal(data, nodata),
    46. origin="upper",
    47. extent=extent,
    48. **kwargs,
    49. )
    50. ax.set_aspect("equal")
    51. ax.set_title(title, fontsize="large")
    52. ax.text(
    53. 0.01, 0.01, "created with pyflwdir", transform=ax.transAxes, fontsize="large"
    54. )
    55. if filename:
    56. plt.savefig(f"{filename}.png")
    57. return ax
    58. # convenience method for vectorizing a raster
    59. def vectorize(data, nodata, transform, crs=crs, name="value"):
    60. feats_gen = features.shapes(
    61. data,
    62. mask=data != nodata,
    63. transform=transform,
    64. connectivity=8,
    65. )
    66. feats = [
    67. {"geometry": geom, "properties": {name: val}} for geom, val in list(feats_gen)
    68. ]
    69. # parse to geopandas for plotting / writing to file
    70. gdf = gpd.GeoDataFrame.from_features(feats, crs=crs)
    71. gdf[name] = gdf[name].astype(data.dtype)
    72. return gdf

    执行代码如下:

    1. import rasterio
    2. import numpy as np
    3. import pyflwdir
    4. from utils import (
    5. quickplot,
    6. colors,
    7. cm,
    8. plt,
    9. ) # data specific quick plot convenience method
    10. # read elevation data of the rhine basin using rasterio
    11. with rasterio.open(r"D:\Raster.tif", "r") as src:
    12. elevtn = src.read(1)
    13. nodata = src.nodata
    14. transform = src.transform
    15. crs = src.crs
    16. extent = np.array(src.bounds)[[0, 2, 1, 3]]
    17. latlon = src.crs.is_geographic
    18. prof = src.profile
    19. ax = quickplot(title="Elevation")
    20. im = ax.imshow(
    21. np.ma.masked_equal(elevtn, -9999),
    22. extent=extent,
    23. cmap="gist_earth_r",
    24. alpha=0.5,
    25. vmin=0,
    26. vmax=1000,
    27. )
    28. fig = plt.gcf()
    29. cax = fig.add_axes([0.8, 0.37, 0.02, 0.12])
    30. fig.colorbar(im, cax=cax, orientation="vertical", extend="max")
    31. cax.set_ylabel("elevation [m+EGM96]")
    32. # plt.savefig('elevation.png', dpi=225, bbox_axis='tight')

    1. #####划分子流域
    2. import geopandas as gpd
    3. import numpy as np
    4. import rasterio
    5. import pyflwdir
    6. # local convenience methods (see utils.py script in notebooks folder)
    7. from utils import vectorize # convenience method to vectorize rasters
    8. from utils import quickplot, colors, cm # data specific quick plot method
    9. # returns FlwDirRaster object
    10. flw = pyflwdir.from_dem(
    11. data=elevtn,
    12. nodata=src.nodata,
    13. transform=transform,
    14. latlon=latlon,
    15. outlets="min",
    16. )
    17. import geopandas as gpd
    18. feats = flw.streams(min_sto=4)
    19. gdf = gpd.GeoDataFrame.from_features(feats, crs=crs)
    20. # create nice colormap of Blues with less white
    21. cmap_streams = colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7)))
    22. gdf_plot_kwds = dict(column="strord", cmap=cmap_streams)
    23. # plot streams with hillshade from elevation data (see utils.py)
    24. ax = quickplot(
    25. gdfs=[(gdf, gdf_plot_kwds)],
    26. title="Streams based steepest gradient algorithm",
    27. filename="flw_streams_steepest_gradient",
    28. )

    1. # update data type and nodata value properties which are different compared to the input elevation grid and write to geotif
    2. prof.update(dtype=d8_data.dtype, nodata=247)
    3. with rasterio.open(r"D:\flwdir.tif", "w", **prof) as src:
    4. src.write(d8_data, 1)
    5. ######################
    6. #Delineation of (sub)basins#
    7. ###########
    8. with rasterio.open(r"D:\flwdir.tif", "r") as src:
    9. flwdir = src.read(1)
    10. crs = src.crs
    11. flw = pyflwdir.from_array(
    12. flwdir,
    13. ftype="d8",
    14. transform=src.transform,
    15. latlon=crs.is_geographic,
    16. cache=True,
    17. )
    18. # define output locations
    19. #x, y = np.array([106.6348,26.8554]), np.array([107.0135,26.8849])
    20. #x, y = np.array([26.8554,106.6348]), np.array([26.8849,107.0135])
    21. x, y = np.array([106.4244,107.0443]), np.array([26.8554,27.0384])
    22. gdf_out = gpd.GeoSeries(gpd.points_from_xy(x, y, crs=4326))
    23. # delineate subbasins
    24. subbasins = flw.basins(xy=(x, y), streams=flw.stream_order() >= 5)
    25. # vectorize subbasins using the vectorize convenience method from utils.py
    26. gdf_bas = vectorize(subbasins.astype(np.int32), 0, flw.transform, name="basin")
    27. gdf_bas.head()
    28. # plot
    29. # key-word arguments passed to GeoDataFrame.plot()
    30. gpd_plot_kwds = dict(
    31. column="basin",
    32. cmap=cm.Set3,
    33. legend=True,
    34. categorical=True,
    35. legend_kwds=dict(title="Basin ID [-]"),
    36. alpha=0.5,
    37. edgecolor="black",
    38. linewidth=0.8,
    39. )
    40. points = (gdf_out, dict(color="red", markersize=20))
    41. bas = (gdf_bas, gpd_plot_kwds)
    42. # plot using quickplot convenience method from utils.py
    43. ax = quickplot([bas, points], title="Basins from point outlets", filename="flw_basins")

    1. # calculate subbasins with a minimum stream order 7 and its outlets
    2. subbas, idxs_out = flw.subbasins_streamorder(min_sto=7, mask=None)
    3. # transfrom map and point locations to GeoDataFrames
    4. gdf_subbas = vectorize(subbas.astype(np.int32), 0, flw.transform, name="basin")
    5. gdf_out = gpd.GeoSeries(gpd.points_from_xy(*flw.xy(idxs_out), crs=4326))
    6. # plot
    7. gpd_plot_kwds = dict(
    8. column="basin", cmap=cm.Set3, edgecolor="black", alpha=0.6, linewidth=0.5
    9. )
    10. bas = (gdf_subbas, gpd_plot_kwds)
    11. points = (gdf_out, dict(color="k", markersize=20))
    12. title = "Subbasins based on a minimum stream order"
    13. ax = quickplot([bas, points], title=title, filename="flw_subbasins")

    1. # get the first level nine pfafstetter basins
    2. pfafbas1, idxs_out = flw.subbasins_pfafstetter(depth=1)
    3. # vectorize raster to obtain polygons
    4. gdf_pfaf1 = vectorize(pfafbas1.astype(np.int32), 0, flw.transform, name="pfaf")
    5. gdf_out = gpd.GeoSeries(gpd.points_from_xy(*flw.xy(idxs_out), crs=4326))
    6. gdf_pfaf1.head()
    7. # plot
    8. gpd_plot_kwds = dict(
    9. column="pfaf",
    10. cmap=cm.Set3_r,
    11. legend=True,
    12. categorical=True,
    13. legend_kwds=dict(title="Pfafstetter \nlevel 1 index [-]", ncol=3),
    14. alpha=0.6,
    15. edgecolor="black",
    16. linewidth=0.4,
    17. )
    18. points = (gdf_out, dict(color="k", markersize=20))
    19. bas = (gdf_pfaf1, gpd_plot_kwds)
    20. title = "Subbasins based on pfafstetter coding (level=1)"
    21. ax = quickplot([bas, points], title=title, filename="flw_pfafbas1")
    1. # lets create a second pfafstetter layer with a minimum subbasin area of 5000 km2
    2. pfafbas2, idxs_out = flw.subbasins_pfafstetter(depth=2, upa_min=5000)
    3. gdf_pfaf2 = vectorize(pfafbas2.astype(np.int32), 0, flw.transform, name="pfaf2")
    4. gdf_out = gpd.GeoSeries(gpd.points_from_xy(*flw.xy(idxs_out), crs=4326))
    5. gdf_pfaf2["pfaf"] = gdf_pfaf2["pfaf2"] // 10
    6. gdf_pfaf2.head()
    7. # plot
    8. bas = (gdf_pfaf2, gpd_plot_kwds)
    9. points = (gdf_out, dict(color="k", markersize=20))
    10. title = "Subbasins based on pfafstetter coding (level=2)"
    11. ax = quickplot([bas, points], title=title, filename="flw_pfafbas2")

    1. # calculate subbasins with a minimum stream order 7 and its outlets
    2. min_area = 2000
    3. subbas, idxs_out = flw.subbasins_area(min_area)
    4. # transfrom map and point locations to GeoDataFrames
    5. gdf_subbas = vectorize(subbas.astype(np.int32), 0, flw.transform, name="basin")
    6. # randomize index for visualization
    7. basids = gdf_subbas["basin"].values
    8. gdf_subbas["color"] = np.random.choice(basids, size=basids.size, replace=False)
    9. # plot
    10. gpd_plot_kwds = dict(
    11. column="color", cmap=cm.Set3, edgecolor="black", alpha=0.6, linewidth=0.5
    12. )
    13. bas = (gdf_subbas, gpd_plot_kwds)
    14. title = f"Subbasins based on a minimum area of {min_area} km2"
    15. ax = quickplot([bas], title=title, filename="flw_subbasins_area")

  • 相关阅读:
    单词相似性查询易语言代码
    Python实现贝叶斯岭回归模型(BayesianRidge算法)并使用K折交叉验证进行模型评估项目实战
    基于RetinaFace的口罩人脸检测算法
    Web渗透之域名(子域名)收集方法
    TaiShan 200服务器安装Ubuntu 18.04
    函数对象的介绍及使用
    DSPE-PEG-PDP,DSPE-PEG-OPSS,磷脂-聚乙二醇-巯基吡啶可减少肽的免疫原性
    大佬的职场经验
    问chatgpt最近生活的困难
    报表工具怎么选?JAVA开源工具那么好用,为什么大家还花钱买商用
  • 原文地址:https://blog.csdn.net/lg8883573/article/details/134500841