我正在按照 Planet 的文档对 Planet Scope 4 波段图像执行 NDVI 计算
以下代码块是我写的：
从工作目录中的原始图像中提取波段数据
import rasterio import numpy
image_file = "20170430_194027_0c82_3B_AnalyticMS"
with rasterio.open(image_file) as src: band_red = src.read(3)
with rasterio.open(image_file) as src: band_nir = src.read(4)
from xml.dom import minidom
xmldoc = minidom.parse("20170430_194027_0c82_3B_AnalyticMS_metadata") nodes = xmldoc.getElementsByTagName("ps:bandSpecificMetadata")
从目录中的元数据文件中提取 TOA 校正系数
TOA_coeffs = {} for node in nodes: bn = node.getElementsByTagName("ps:bandNumber")[0].firstChild.data if bn in ['1', '2', '3', '4']:
i = int(bn)
value = node.getElementsByTagName("ps:ReflectanceCoefficient")[0].firstChild.data
TOA_coeffs[1] = float(value)
计算 NDVI 并保存文件
band_red = band_red * TOA_coeffs[3] band_nir = band_nir * TOA_coeffs[4]
numpy.seterr(divide = 'ignore', invalid = 'ignore')
NDVI = (band_nir.astype(float) - band_red.astype(float))/(band_nir + band_red) numpy.nanmin(NDVI), numpy.nanmax(NDVI)
kwargs = src.meta kwargs.update(dtype=rasterio.float32,
count = 1)
with rasterio.open('ndvi.tif', 'W', **kwargs) as dst: dst.write_band(1, NDVI.astype(rasterio.float32))
添加符号系统并绘制颜色条
import matplotlib.pyplot as plt import matplotlib.colors as colors
class MidpointNormalize(colors.Normalize): def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
self.midpoint = midpoint
colors.Normalize.__init__(self, vmin, vmax, clip)
def __call__(self, value, clip=None):
x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
return numpy.ma.masked_array(numpy.interp(value, x, y), >numpy.isnan(value))
min = numpy.nanmin(NDVI) min = numpy.nanmax(NDVI) mid = 0.1
fig = plt.figure(figsize= (20,10)) ax = fig.add_subplot(111)
cmap = plt.cm.RdYlGn
cax = ax.imshow(NDVI, cmap=cmap, clim=(min,max),
>norm=MidpointNormalize(midpoint=mid, vmin=min, vmax=max))
ax.axis('off') ax.set_title('NDVI_test', fontsize= 18, fontweight='bold')
cbar = fig.colorbar(cax, orientation= 'horizontal', shrink=0.65)
fig.savefig("output/NDVI_test.png", dpi=200, bbox_inches='tight',
>pad_inches=0.7)
plt.show()
绘制 NDVI 像素值分布的直方图
fig2 = plt.figure(figsize=(10,10)) ax = fig2.add_subplot(111)
plt.title("NDVI Histogram", fontsize=18, fontweight='bold') plt.xlabel("NDVI values", fontsize=14) plt.ylabel("# pixels", fontsize=14)
x = NDVI[~numpy.isnan(NDVI)] numBins = 20 ax.hist(x,numBins,color='green',alpha=0.8)
fig2.savefig("output/ndvi-histogram.png", dpi=200, bbox_inches='tight', >pad_inches=0.7)
plt.show()
唉，脚本的执行在代码的开头被缩短了：
File "C:/Users/David/Desktop/ArcGIS files/Planet Labs/2017.6_Luis_Bedin_Bolivia/planet_order_58311/20170430_194027_0c82/TOA_correction_NDVI.py", line 8, in <module>
import rasterio
ModuleNotFoundError: No module named 'rasterio'
所以我决定安装光栅，这应该可以解决问题：
C:\Users\David\Desktop\ArcGIS files\Planet Labs\2017.6_Luis_Bedin_Bolivia\planet_order_58311\20170430_194027_0c82>pip install rasterio
Collecting rasterio
Using cached rasterio-0.36.0.tar.gz
Requirement already satisfied: affine in c:\users\david\anaconda3\lib\site-packages (from rasterio)
Requirement already satisfied: cligj in c:\users\david\anaconda3\lib\site-packages (from rasterio)
Requirement already satisfied: numpy in c:\users\david\anaconda3\lib\site-packages (from rasterio)
Requirement already satisfied: snuggs in c:\users\david\anaconda3\lib\site-packages (from rasterio)
Requirement already satisfied: click-plugins in c:\users\david\anaconda3\lib\site-packages (from rasterio)
我从中解释的是已经安装了 rasterio。如果 Python 控制台告诉我没有名为 rasterio 的模块，这怎么可能呢？控制台的输出还显示需要 Microsoft Visual C++。经过进一步研究，我找到了这个用户的解决方案。我试过了，但控制台也告诉我已经安装了 rasterio：
(envpythonfs) C:\Users\David\Desktop\ArcGIS files\Planet Labs\2017.6_Luis_Bedin_Bolivia\planet_order_58311\20170430_194027_0c82>conda install rasterio gdal
Fetching package metadata .............
Solving package specifications: .
# All requested packages already installed.
# packages in environment at C:\Users\David\Anaconda3\envs\envpythonfs:
#
我在 Windows 10 64 位机器上使用 Spyder 3.1.2 和 Python 3.6 创建脚本。