1

python - 如何在Python中沿给定的shapefile行从栅格中提取值的配置文件？

我正在努力寻找一种从栅格（geotiff）中提取值轮廓（例如地形轮廓）的方法。库Rasterio 有一种方法可以从基于多边形的栅格中剪辑/提取值，但我找不到线 shapefile 的等效方法。

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

换句话说，我正在寻找与 QGIS 中的 Terrain Profile 工具提供的 Python 等效项。

谢谢

4

2 回答 2

1

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

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

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

如果您只想获取该行下的原始像素值，您应该能够为每一行使用掩码rasterio或直接光栅化。您可能想all_touched=True在线条的情况下使用。

于 2020-06-10T14:18:09.747 回答
0

我遇到了类似的问题，并找到了适合我的解决方案。该解决方案用于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()
于 2021-12-20T15:57:27.597 回答