本文介绍了如何使用Python编程提取栅格图像（如Landsat影像生成的NDVI图像）中，特定多边形区域内部和外部的NDVI均值。通过rasterio和fiona库，可以轻松实现对栅格数据的裁剪、掩膜以及统计计算，为遥感数据分析提供了一种高效且可重复的方法。
准备工作
在开始之前，需要确保安装了以下Python库：
rasterio: 用于读取和写入栅格数据。fiona: 用于读取矢量数据（例如，shapefile）。numpy: 用于数值计算，例如计算均值。
可以使用pip命令安装这些库：
pip install rasterio fiona numpy
核心步骤
读取矢量数据 (Shapefile)
立即学习“Python免费学习笔记（深入）”；
首先，需要使用fiona库读取包含多边形信息的shapefile。以下代码展示了如何读取shapefile并提取几何信息：
import fionashapefile_path = "path/to/your/shapefile.shp" # 替换为你的shapefile路径with fiona.open(shapefile_path, "r") as sf: shapes = [feature["geometry"] for feature in sf]
这段代码打开指定的shapefile，并提取每个要素的几何信息，存储在shapes列表中。shapes列表将包含代表多边形的几何对象。
读取栅格数据 (TIFF)
接下来，使用rasterio库读取栅格数据（例如，NDVI图像）。以下代码展示了如何读取TIFF图像：
import rasterioraster_path = "path/to/your/ndvi.tif" # 替换为你的NDVI图像路径with rasterio.open(raster_path) as src: raster_image = src.read(1) # 读取第一个波段，假设NDVI数据存储在第一个波段 raster_meta = src.meta.copy() # 复制元数据，后续写入使用
这段代码打开指定的TIFF图像，并将图像数据读取到raster_image变量中。同时，复制图像的元数据，以便后续写入处理后的数据。
使用多边形裁剪栅格数据
使用rasterio.mask.mask函数，可以根据多边形裁剪栅格数据。以下代码展示了如何裁剪栅格数据，并计算裁剪区域内的NDVI均值：
import rasterio.maskimport numpy as npwith rasterio.open(raster_path) as src: out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True) out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致# 计算裁剪区域内的NDVI均值，忽略无效值NDVI_mean = np.nanmean(out_image)print(f"多边形内部NDVI均值: {NDVI_mean}")
这段代码使用shapes列表中的多边形裁剪栅格图像。crop=True参数表示裁剪到最小范围。np.nanmean用于计算均值，并忽略NaN值。
提取多边形外部的NDVI值
可以使用rasterio.mask.mask函数，通过设置invert=True参数来提取多边形外部的像素值。以下代码展示了如何提取多边形外部的NDVI值：
with rasterio.open(raster_path) as src: out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True, invert=True) out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致# 计算裁剪区域外的NDVI均值，忽略无效值NDVI_mean_outside = np.nanmean(out_image)print(f"多边形外部NDVI均值: {NDVI_mean_outside}")
通过设置invert=True，rasterio.mask.mask函数将返回多边形外部的像素值。
分别处理多个多边形
如果需要分别处理多个多边形，可以循环遍历shapes列表，并对每个多边形执行裁剪和均值计算操作。
for i, shape in enumerate(shapes): with rasterio.open(raster_path) as src: out_image, out_transform = rasterio.mask.mask(src, [shape], crop=True) out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致 NDVI_mean = np.nanmean(out_image) print(f"多边形 {i+1} 内部NDVI均值: {NDVI_mean}") with rasterio.open(raster_path) as src: out_image, out_transform = rasterio.mask.mask(src, [shape], crop=True, invert=True) out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致 NDVI_mean_outside = np.nanmean(out_image) print(f"多边形 {i+1} 外部NDVI均值: {NDVI_mean_outside}")
这段代码循环遍历shapes列表，并对每个多边形分别计算内部和外部的NDVI均值。
完整示例代码
import rasterioimport fionaimport rasterio.maskimport numpy as npshapefile_path = "path/to/your/shapefile.shp" # 替换为你的shapefile路径raster_path = "path/to/your/ndvi.tif" # 替换为你的NDVI图像路径# 读取矢量数据with fiona.open(shapefile_path, "r") as sf: shapes = [feature["geometry"] for feature in sf]# 循环处理每个多边形for i, shape in enumerate(shapes): # 多边形内部 with rasterio.open(raster_path) as src: out_image, out_transform = rasterio.mask.mask(src, [shape], crop=True) out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致 NDVI_mean = np.nanmean(out_image) print(f"多边形 {i+1} 内部NDVI均值: {NDVI_mean}") # 多边形外部 with rasterio.open(raster_path) as src: out_image, out_transform = rasterio.mask.mask(src, [shape], crop=True, invert=True) out_image = out_image.astype(src.dtypes[0]) # 确保数据类型一致 NDVI_mean_outside = np.nanmean(out_image) print(f"多边形 {i+1} 外部NDVI均值: {NDVI_mean_outside}")
注意事项
确保shapefile和栅格图像的坐标系一致。如果坐标系不一致，需要进行投影转换。rasterio.mask.mask函数返回的out_image是一个三维数组，第一个维度的大小为1，需要使用out_image[0]来访问二维图像数据。在计算均值时，使用np.nanmean可以忽略无效值（例如，NaN值）。根据实际需求，可以修改代码以提取其他统计信息，例如最大值、最小值、标准差等。需要确保栅格数据类型与计算时的数据类型一致，如果数据类型不一致，可能会导致计算结果错误，使用 out_image = out_image.astype(src.dtypes[0]) 可以保证数据类型一致。
总结
本文介绍了如何使用Python的rasterio和fiona库提取栅格数据中多边形内外NDVI均值。通过这些方法，可以高效地进行遥感数据分析，并提取有用的信息。 掌握这些技术，可以为更复杂的遥感应用打下坚实的基础。
以上就是利用Python提取栅格数据中多边形内外NDVI均值的详细内容，更多请关注创想鸟其它相关文章！
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。
如发现本站有涉嫌抄袭侵权/违法违规的内容， 请发送邮件至 chuangxiangniao@163.com 举报，一经查实，本站将立刻删除。
发布者：程序猿，转转请注明出处：https://www.chuangxiangniao.com/p/1367653.html
微信扫一扫 
支付宝扫一扫 
