精准农业中遥感技术应用（四）- 作物长势指数计算演示

本文链接：https://blog.csdn.net/CobyBryantGiGi/article/details/142439180

上篇文章介绍一些重要的分析作物长势的指数，本章使用 EarthPy包分析某地理区间内作物长势的情况。

EarthPy介绍

EarthPy 是一个开源的 Python 库，旨在简化地球科学数据的处理、分析和可视化。它建立在其他强大的开源库之上，如 GeoPandas、Rasterio、Matplotlib 和 Pandas，为地球科学领域的研究人员、教育工作者和学生提供了一个方便的工具集。

功能列表

数据处理与可视化：
- EarthPy 提供了一套工具，可以轻松加载和处理栅格和矢量数据。
- 它支持数据标准化、阴影计算和空间范围计算等任务，这些都是环境科学研究中常见的需求。
遥感数据分析：
- 通过 EarthPy，用户可以方便地查看遥感图像中单个波段的直方图，这对于数据校准和质量控制非常重要。
- 它还允许用户创建底图、绘制图像和渲染多波段光谱遥感图像，这些功能对于数据的可视化和解释非常有用。
教育与教学：
- EarthPy 的 io 模块特别有助于教育领域，它简化了教学数据的下载和共享，使得课堂活动更加高效。
- 它支持工作流程的可重复性，这对于教学和学习都是至关重要的。
空间分析：
- EarthPy 支持空间数据的输入/输出、掩膜处理、绘图和空间分析，这些功能对于地理空间数据的操作至关重要。
- 它提供了一个高效的平台，使得研究人员能够快速进行空间数据的探索和清理。
- 通过集成 GeoPandas 和 Rasterio，EarthPy 允许用户在同一个包中处理矢量和栅格数据，这大大提高了工作效率。

优势、特点

EarthPy在遥感数据分析和应用上的优势和特点包括：

简化数据处理：提供易于使用的函数来处理栅格和矢量数据，简化了数据加载、预处理和分析的步骤。
直方图和图形探索：允许用户查看遥感图像中单个波段的直方图，帮助探索数据校准和质量问题。
多波段光谱渲染：支持创建具有独特符号的底图和绘制带有颜色条的图像，便于多波段光谱遥感图像的可视化。
空间分析工具：集成了空间分析功能，如掩膜处理、绘图和空间分析，提高了地理空间数据操作的效率。
教育和研究支持：特别适合教育和研究环境，因为它提供了丰富的文档和教程，有助于学习和应用遥感数据分析。
开源和社区驱动：作为一个开源项目，EarthPy受益于活跃的社区贡献，不断更新和改进，以满足用户的需求。

作物长势指数演算示例

统一植被指数（NDVI）计算

import os
from glob import glob
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep

# Get data and set your home working directory
data = et.data.get_data("vignette-landsat")

os.chdir(os.path.join(et.io.HOME, '/Users/****/T-Code/earthpy-data/earth-analytics', 'data'))

# Stack the Landsat 8 bands
# This creates a numpy array with each "layer" representing a single band
landsat_path = glob( "/Users/****/earthpy-data/earth-analytics/data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*_crop.tif")
landsat_path.sort()

# print(f"{landsat_path}");
arr_st, meta = es.stack(landsat_path, nodata=-9999)

# Landsat 8 red band is band 4 at [3]
# Landsat 8 near-infrared band is band 5 at [4]
ndvi = es.normalized_diff(arr_st[4], arr_st[3])
titles = ["Normalized Difference Vegetation Index (NDVI)"]
# Turn off bytescale scaling due to float values for NDVI
ep.plot_bands(ndvi, cmap="RdYlGn", cols=1, title=titles, vmin=-1, vmax=1)

代码说明：

下载Landsat的卫星遥感示例数据，并设置本地的操作环境，然后读取不同波段的多光谱数据，使用近红外光和红光波段，计算归一化统一植被指数，并绘制图片展示。注意，运行过程中可能因为工程目录设置，导致找不到文件的情况，可以设置多光谱存放目录的绝对路径。

运行结果：

统一植被指数的分类

ndvi_class_bins = [-np.inf, 0, 0.15, 0.23, 0.6, np.inf]

ndvi_landsat_class = np.digitize(ndvi, ndvi_class_bins)

# Apply the nodata mask to the newly classified NDVI data
ndvi_landsat_class = np.ma.masked_where(np.ma.getmask(ndvi), ndvi_landsat_class)
np.unique(ndvi_landsat_class)
nbr_colors = ["gray", "y", "yellowgreen", "g", "darkgreen"]
nbr_cmap = ListedColormap(nbr_colors)

# Define class names
ndvi_cat_names = [
    "No Vegetation",
    "Bare Area",
    "Low Vegetation",
    "Moderate Vegetation",
    "High Vegetation",
]

# Get list of classes
classes = np.unique(ndvi_landsat_class)
classes = classes.tolist()

# The mask returns a value of none in the classes. remove that
classes = classes[0:5]

# Plot the data
fig, ax = plt.subplots(figsize=(10, 10))

im = ax.imshow(ndvi_landsat_class, cmap=nbr_cmap)

epp.draw_legend(im_ax=im, classes=classes, titles=ndvi_cat_names)

ax.set_title(
    "Landsat 8 - Normalized Difference Vegetation Index (NDVI) Classes",
    fontsize=14,
)
ax.set_axis_off()

# Auto adjust subplot to fit figure size
plt.tight_layout()

运行结果：