为了账号安全,请及时绑定邮箱和手机立即绑定

如何沿给定线从栅格中提取价值概况?

如何沿给定线从栅格中提取价值概况?

沧海一幻觉 2023-01-04 16:20:35

如何在 Python 中沿着给定的 shapefile 线从栅格中提取值的配置文件?

我正在努力寻找一种方法来从栅格 (geotiff) 中提取值的配置文件(例如地形配置文件)。库Rasterio 有一种方法可以根据多边形从栅格中裁剪/提取值,但我找不到线形文件的等效方法。

scipy有一个基本方法,但它并不像 rasterio 可以提供的基于更高级别工具箱的方法那样固有地保存地理信息。

换句话说,我正在寻找与 QGIS 中的 Terrain Profile 工具相同的 Python 工具。


    查看完整描述

    2 回答

    ?
    芜湖不芜

    TA贡献1544条经验 获得超7个赞

    这与提取多边形有点不同,因为您想要按照触摸的顺序对线触摸的每个像素进行采样(多边形方法不关心像素顺序)。

    看起来可以改用这种方法rasterio代替使用。geopandas给定使用或fiona作为对象从 shapefile 中读取的一条线shapely,您可以使用端点导出一个新的等距投影,您可以dst_crsWarpedVRT中使用该投影并从中读取像素值。看起来你需要根据你想要采样的像素数来计算你的线的长度,这是WarpedVRT.

    如果您的线不是端点之间的近似直线,则可能需要进一步调整此方法。

    如果您只想获取线下的原始像素值,您应该能够为每条线使用遮罩rasterio或直接栅格化。对于线条,您可能想要使用all_touched=True


    查看完整回答
    反对 回复 2023-01-04
    ?
    哈士奇WWW

    TA贡献1547条经验 获得超6个赞

    我遇到了类似的问题,并找到了适合我的解决方案。该解决方案用于shapely对一条/多条线上的点进行采样,然后从 GeoTiff 访问相应的值,因此提取的配置文件遵循线的方向。这是我最终得到的方法:


    def extract_along_line(xarr, line, n_samples=256):

        profile = []


        for i in range(n_samples):

            # get next point on the line

            point = line.interpolate(i / n_samples - 1., normalized=True)

            # access the nearest pixel in the xarray

            value = xarr.sel(x=point.x, y=point.y, method="nearest").data

            profile.append(value)

            

        return profile

    这是一个使用数据库中的数据的工作示例,copernicus-dem该线是接收到的图块的对角线:


    import rioxarray

    import shapely.geometry

    import matplotlib.pyplot as plt


    sample_tif = ('https://elevationeuwest.blob.core.windows.net/copernicus-dem/'

                  'COP30_hh/Copernicus_DSM_COG_10_N35_00_E138_00_DEM.tif')


    # Load xarray

    tile = rioxarray.open_rasterio(sample_tif).squeeze()

    # create a line (here its the diagonal of tile)

    line = shapely.geometry.MultiLineString([[

                [tile.x[-1],tile.y[-1]],

                [tile.x[0], tile.y[0]]]])


    # use the method from above to extract the profile

    profile = extract_along_line(tile, line)

    plt.plot(profile)

    plt.show()


    查看完整回答
    反对 回复 2023-01-04

    添加回答

    举报

    0/150
    提交
    取消
    意见反馈 帮助中心 APP下载
    官方微信