如何计算 GRASS GIS 中未嵌套的流域？答案

【问题标题】：How to to calculate unnested watersheds in GRASS GIS?如何计算 GRASS GIS 中未嵌套的流域？
【发布时间】：2023-02-15 18:41:53
【问题描述】：

在 Python 中运行 GRASS GIS 模块 r.accumulate 时，我遇到了一些问题。我使用该模块计算 7000 多个测量点的子流域。不幸的是，算法的输出是嵌套的。所以所有的子流域都相互重叠。运行 r.accumulate 子流域模块对于一个或多个点大约需要 2 分钟，我假设瓶颈是加载方向栅格。

我想知道 GRASS GIS 中是否有可用的未嵌套变体，如果没有，如何克服每次调用模块 accumulate 时加载方向栅格的瓶颈。下面是我到目前为止尝试过的代码 sn-p（导致嵌套变体）：

locations = VectorTopo('locations',mapset='PERMANENT')
    locations.open('r')
    points=[]
    for i in range(len(locations)):
        points.append(locations.read(i+1).coords())
    for j in range(0,len(points),255):
        output = "watershed_batch_{}@Watersheds".format(j)
        gs.run_command("r.accumulate", direction='direction@PERMANENT', subwatershed=output,overwrite=True, flags = "r", coordinates = points[j:j+255])
        gs.run_command('r.stats', flags="ac", input=output, output="stat_batch_{}.csv".format(j),overwrite=True)

任何想法或想法都非常受欢迎。

【问题讨论】：

标签： python watershed grass

【解决方案1】：

我已经回复了你的邮件，但现在我看到了你的 Python 代码并且更好地理解了你的“重叠”问题。在这种情况下，您不想一次为单个插座点供电。你可以跑

r.accumulate direction=direction@PERMANENT subwatershed=output outlet=locations

r.accumulate 的 outlet 选项可以处理多个 outlet 并将生成不重叠的子流域。

【讨论】：

【解决方案2】：

通过电子邮件提供的答案非常有用。为了分享答案，我提供了下面的代码来进行未嵌套的流域子流域计算。一个小评论：我不得不分批输入坐标，因为坐标列表超出了窗口可以处理的最大字符长度。

感谢@Huidae Cho，调用 R.accumulate 来计算子流域和最长流路径现在可以在一个调用中完成，而不是两个单独的调用。

输出是未嵌套的盆地。较大的子流域与较小的子流域分开，而不是被剪裁成较小的流域。这与输出是栅格格式这一事实有关，其中每个像元只能代表一个盆地。

gs.run_command('g.mapset',mapset='Watersheds')
gs.run_command('g.region', rast='direction@PERMANENT')
StationIds = list(gs.vector.vector_db_select('locations_snapped_new', columns = 'StationId')["values"].values())
XY = list(gs.vector.vector_db_select('locations_snapped_new', columns = 'x_real,y_real')["values"].values())

for j in range(0,len(XY),255):
    output_ws = "watershed_batch_{}@Watersheds".format(j)
    output_lfp = "lfp_batch_{}@Watersheds".format(j)
    output_lfp_unique = "lfp_unique_batch_{}@Watersheds".format(j)
    gs.run_command("r.accumulate", direction='direction@PERMANENT', subwatershed=output_ws, flags = "ar", coordinates = XY[j:j+255],lfp=output_lfp, id=StationIds[j:j+255], id_column="id",overwrite=True)
    gs.run_command("r.to.vect", input=output_ws, output=output_ws, type="area", overwrite=True)
    gs.run_command("v.extract", input=output_lfp, where="1 order by id", output=output_lfp_unique,overwrite=True)

为了导出独特的流域，我使用了以下代码。我必须将 longest_flow_path 转换为指向一些 longest_flow_paths 与其旁边的分水岭的角边界相交。因此，一些最长的水流路径并未完全位于次流域内。请参见下图，其中红线（最长的水流路径）接触子流域边界： enter image description here

gs.run_command('g.mapset',mapset='Watersheds')
lfps= gs.list_grouped('vect', pattern='lfp_unique_*')['Watersheds']
ws= gs.list_grouped('vect', pattern='watershed_batch*')['Watersheds']
files=np.stack((lfps,ws)).T
#print(files)
for file in files:
    print(file)
    ids = list(gs.vector.vector_db_select(file[0],columns="id")["values"].values())
    for idx in ids:
        idx=int(idx[0])
        expr = f'id="{idx}"'
        gs.run_command('v.extract',input=file[0], where=expr, output="tmp_lfp",overwrite=True)
        gs.run_command("v.to.points", input="tmp_lfp", output="tmp_lfp_points", use="vertex", overwrite=True)
        gs.run_command('v.select', ainput= file[1], binput = "tmp_lfp_points", output="tmp_subwatersheds", overwrite=True)
        gs.run_command('v.db.update', map = "tmp_subwatersheds",col= "value", value=idx)
        gs.run_command('g.mapset',mapset='vector_out')
        gs.run_command('v.dissolve',input= "tmp_subwatersheds@Watersheds", output="subwatersheds_{}".format(idx),col="value",overwrite=True)
        gs.run_command('g.mapset',mapset='Watersheds')
        gs.run_command("g.remove", flags="f", type="vector",name="tmp_lfp,tmp_subwatersheds")

我最终得到了每个子流域的向量

【讨论】：