python 提取 栅格tif |
您所在的位置：网站首页 › python读取tif栅格数据 › python 提取 栅格tif |
python 提取 栅格tif
|
在进行地理数据处理的时候，我们经常需要用点矢量来提取栅格数据中的值，在Arcgis或者R里面，这个很好实现，用extract工具或者函数就行，但是我看了下貌似Python里面没有现成的。那就自己写一个吧。 本文示例用一个点矢量的shp文件来批量提取一个栅格数据上各点处的值。 公众号《GIS与Climate》，欢迎关注和交流。 读取数据分别用rasterio和geopandas包读取栅格和矢量数据： rs = rasterio.open('../input/rs_min_1km.tif')shp = gdp.read_file('../input/shp/test_points_proj.shp')
image
可视化一下： fig, ax = plt.subplots()shp.plot(ax=ax,color='orangered')show(rs,ax=ax)
叠加显示
shp文件通过geopandas读取进来之后是一个geodataframe格式的数据，不要多想，直接把它看做是一个dataframe格式的就行（类似于R语言中的sf包），这样子pandas dataframe上可以用的方法我们也可以直接用在geodataframe上，那么批量处理就方便了。 从栅格提取值的本质要知道，栅格数据在本质上就是一个array，如果想要知道某个点的值，其实就是获取这个array中某一个元素的值。如果是一个简单的numpy array，那么知道行列号就可以直接通过索引来获取值了，但是对于地理数据，我们一般是根据地理坐标来提取对应位置的栅格值，而不是行列号。 知道了本质就好办了，把对应的地理坐标变换到矩阵对应的行列号就可以了，本质上这是两个坐标空间的变换，也就是仿射变换。大佬们已经写好了这个功能，直接用rasterio包中的index方法即可，具体见参考链接【1】。 自定义提取函数然后我们根据上面的原理写一个函数： def ExtractPointValue(geometry,rs): ''' geometry: geodataframe的geometry列 rs: 栅格数据 ''' x = geometry.xy[0][0] y = geometry.xy[1][0] row, col = rs.index(x,y) value = rs.read(1)[row,col] return value 批量提取值有了函数之后，我们用apply函数直接把上述函数作用于geodataframe的geometry列即可： shp['value1'] = shp['geometry'].apply(ExtractPointValue,rs=rs) shp.head()
可以看到已经成功提取了各个点处的值。 总结 通过点来提取栅格处的值本质就是进行矩阵元素的索引； 巧妙利用pandas的apply函数进行批量处理； 注意矢量和栅格的坐标系要一致，不一致的需要先统一。 小思考 如果点超级多怎么办？ 栅格运算永远比矢栅叠加运算快，先把点栅格化再获取索引就更快了 如何提取矢量多边形的值？ 下篇文章。 如何批量提取栅格时间序列的值？ 上面的函数改一下，read函数直接读取多波段栅格就行（不加1），或者加个循环； 不是多波段的图像，再用apply函数，换个参数就行；公众号《GIS与Climate》，欢迎关注和交流。 
参考
【1】https://rasterio.readthedocs.io/en/latest/quickstart.html?highlight=index#spatial-indexing 【2】https://hatarilabs.com/ih-en/extract-point-value-from-a-raster-file-with-python-geopandas-and-rasterio-tutorial 【3】https://gis.stackexchange.com/questions/260304/extract-raster-values-within-shapefile-with-pygeoprocessing-or-gdal 【4】https://www.earthdatascience.org/courses/use-data-open-source-python/spatial-data-applications/lidar-remote-sensing-uncertainty/extract-data-from-raster/ |
【本文地址】
今日新闻 |
推荐新闻 |