前言
MODIS 是一种卫星遥感仪器，每天以 250-500 米的分辨率在全球范围内收集数据。 了解如何在 Python 中导入、清理和绘制 MODIS 数据。
1、MODIS 影像简介
中分辨率成像光谱仪 (MODIS) 是一种基于卫星的仪器，可连续收集地球表面的数据。 目前，MODIS 拥有公开可用遥感数据的最佳时间分辨率，每 24 小时覆盖整个地球。
MODIS 收集 36 个光谱波段的数据； 但是，在本博客中，您只会使用前 7 个波段。
（1）MODIS 地表反射率（MOD09GA 产品）
这里有许多不同的 MODIS 数据产品。 这些是经过处理以用于科学的数据集。 在本课程中，我们使用的是 MOD09GA 产品，这是一种包含 MODIS 前 7 个波段的反射率产品。
地表反射率值的正常范围是 0 到 1，其中 1 是最亮的值，0 是最暗的值。 表面反射率是地球表面光谱反射率的量度，就像在地面上测量一样。 你可以把它想象成你的眼睛会看到什么，当然除了你的眼睛看不到电磁波谱可见部分之外的光。
MODIS 提供了许多标准化产品，包括您将在本课程中使用的表面反射率 MOD09GA 产品。 MOD09GA 产品在上表中列出的 7 个光谱带中以 500m 的空间分辨率提供表面反射率。
根据创建 MOD09 产品的 Land Surface Reflectance Science Computing Facility，这些产品是每个波段的表面光谱反射率的估计值，因为它会在地面上测量，就好像没有大气散射或吸收一样。 它校正了大气气体、气溶胶和薄卷云的影响。
（2）MOD09GA 产品的波段元数据
要更好地了解 MODIS 数据，请查看第 14 页 MODIS 用户指南中 MOD09GA 产品的详细表格。
部分表格如下：
使用 MOD09GA 产品的此表，回答以下问题：
我们数据的有效值范围是多少？
什么没有数据值？
与我们的数据相关的比例因子是多少？
（3）识别用于 NBR 计算的 MODIS 波段
本博客上，您将使用 MODIS 数据计算 NBR。 然而，即使您可以使用许多不同的遥感产品计算相同的植被指数，请记住每个遥感数据的波段是不同的。
查看上表，其中显示了 MODIS 传感器的波段范围。 回想一下，NBR 指数适用于 NIR 波段在 760 – 900 nm 和 SWIR 波段在 2080 – 2350 nm 之间的任何多光谱传感器。
2、打开 MODIS 图像
您将学习如何使用科罗拉多州冷泉火灾研究区的火灾前 MODIS 影像打开 MODIS 数据。
在开始之前，导入以下包并确保设置了工作目录。
from glob import glob
import os
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
import geopandas as gpd
import rioxarray as rxr
import xarray as xr
from rasterio.plot import plotting_extent
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
from shapely.geometry import box
data = et.data.get_data('cold-springs-fire')
os.chdir(os.path.join(et.io.HOME, 'earth-analytics', 'data'))
Downloading from https://ndownloader.figshare.com/files/10960109
Extracted output to /root/earth-analytics/data/cold-springs-fire/.在之前的课程中，您已经使用 glob(“*keyword*.tif”) 创建了所有文件的列表，这些文件同时：
包含由星号表示的特定关键字（例如 *band*）和
包含扩展名 .tif。
首先使用 glob 使用关键字和扩展名 *sur_refl_b*.tif 创建 MODIS 表面反射栅格列表。
# Create list of MODIS rasters for surface reflectance
modis_bands_pre_list = glob(os.path.join("cold-springs-fire",
"modis",
"reflectance",
"07_july_2016",
"crop",
"*_sur_refl_b*.tif"))
# Sort the list of bands
modis_bands_pre_list.sort()
modis_bands_pre_list
['cold-springs-fire/modis/reflectance/07_july_2016/crop/MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b01_1.tif',
'cold-springs-fire/modis/reflectance/07_july_2016/crop/MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b02_1.tif',
'cold-springs-fire/modis/reflectance/07_july_2016/crop/MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b03_1.tif',
'cold-springs-fire/modis/reflectance/07_july_2016/crop/MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b04_1.tif',
'cold-springs-fire/modis/reflectance/07_july_2016/crop/MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b05_1.tif',
'cold-springs-fire/modis/reflectance/07_july_2016/crop/MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b06_1.tif',
'cold-springs-fire/modis/reflectance/07_july_2016/crop/MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b07def combine_tifs(tif_list):
"""A function that combines a list of tifs in the same CRS
and of the same extent into an xarray object
Parameters
----------
tif_list : list
A list of paths to the tif files that you wish to combine.
Returns
-------
An xarray object with all of the tif files in the listmerged into
a single object.
"""
out_xr = []
for i, tif_path in enumerate(tif_list):
out_xr.append(rxr.open_rasterio(tif_path, masked=True).squeeze())
out_xr[i]["band"] = i+1
return xr.concat(out_xr, dim="band")接下来，使用函数 combine_tifs 从您使用 glob 创建的波段列表创建 rioxarray 对象。 然后，您可以导入 MODIS 波段并创建 RGB 图。
# Open file list with function
modis_bands_pre = combine_tifs(modis_bands_pre_list)# Plot MODIS RGB
ep.plot_rgb(modis_bands_pre.values,
rgb=[0, 3, 2],
title="Surface Reflectance \n MODIS RGB Bands")
plt.show()
使用 RGB 波段的 MODIS 表面反射率，用于预冷泉火灾时间段。
3、查看数据的值
要开始探索数据，您可以计算选定波段的最小值和最大值以查看值的范围。 例如，您可以为 MODIS 堆栈的第一个波段（红色）计算这些值。
# Identify minimum and maximum values of band 1 (red)
print(modis_bands_pre[1].min(), modis_bands_pre[1].max())
<xarray.DataArray ()>
array(-100., dtype=float32)
Coordinates:
band int64 2
spatial_ref int64 0 <xarray.DataArray ()>
array(10039., dtype=float32)
Coordinates:
band int64 2
spatial_ref int64 0正在打开的数据很棒！ 但是，您还需要为作业裁剪数据。 您可以通过两个简单的步骤完成此操作。
重新投影火灾边界，使其与您的 MODIS 数据位于同一 CRS 中
使用火边界和 xarray_name.rio.clip() 函数裁剪 MODIS 数据。 要将此函数用于 GeoPandas 对象，您必须调用 geopandas 对象的几何列。
xarray_name.rio.clip(crop_bound.geometry)但是，此方法会将场景剪裁成矢量形状的确切轮廓。 这可能很有用，但通常您会希望剪切覆盖矢量形状边界的矩形，而不是精确的形状。 为此，您可以使用从 shapely.geometry 导入的 box 函数在要剪切的几何体周围创建一个框。
xarray_name.rio.clip([box(*crop_bound.total_bounds)])您可以在下面看到这两种方法剪辑场景的不同之处！
# Open fire boundary
fire_boundary_path = os.path.join("cold-springs-fire",
"vector_layers",
"fire-boundary-geomac",
"co_cold_springs_20160711_2200_dd83.shp")
fire_boundary = gpd.read_file(fire_boundary_path)
fire_bound_sin = fire_boundary.to_crs(modis_bands_pre.rio.crs)
fire_bound_box = [box(*fire_bound_sin.total_bounds)]
# MODIS Clipped to Geometry
modis_clip_geometry = modis_bands_pre.rio.clip(fire_bound_sin.geometry,
all_touched=True,
from_disk=True)
# MODIS Clipped to Bounds
modis_clip = modis_bands_pre.rio.clip(fire_bound_box,
all_touched=True,
from_disk=True)
extent = plotting_extent(modis_clip[0].values, modis_clip.rio.transform())
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 15))
# Plotting Geometry Clip
ep.plot_rgb(modis_clip.values,
rgb=[0, 3, 2],
ax=ax1,
extent=extent,
title='MODIS Clipped to Bounds')
fire_bound_sin.boundary.plot(ax=ax1)
# Plotting Bounds Clip
ep.plot_rgb(modis_clip_geometry.values,
rgb=[0, 3, 2],
ax=ax2,
extent=extent,
title='MODIS Clipped to Geometry')
fire_bound_sin.boundary.plot(ax=ax2)
plt.show()
出于 NDVI 和 NBR 分析的目的，我们对周围区域以及直接位于几何体内部的像素感兴趣，因此对于此分析，我们将使用边界方法。
要探索裁剪后的数据，请回想一下，您还可以使用 earthpy 中的 plot_bands 函数来创建每个波段的图。
# Create a list of titles
titles = ["Red Band", "Near Infrared (NIR) Band", "Blue/Green Band", "Green Band",
"Near Infrared (NIR) Band", "Mid-infrared Band", "Mid-infrared Band"]
# Plot all bands individually
ep.plot_bands(modis_clip,
cols=3,
title=titles,
figsize=(10, 6))
plt.show()
来自 MODIS 的表面反射率裁剪图像，适用于冷泉前火灾的所有波段。
使用 RGB 波段的 MODIS 裁剪表面反射率，用于预冷泉火灾。
# Create a colors and titles list to use in the histogram
colors = ['r', 'k', 'b', 'g', 'k', 'y', 'y']
titles = ["Red Band", "Near Infrared (NIR) Band", "Blue/Green Band",
"Green Band", "Near Infrared (NIR) Band",
"Mid-infrared Band", "Mid-infrared Band"]
# Plot histogram
ep.hist(modis_clip.values,
colors=colors,
title=titles,
cols=2)
plt.show()
来自 MODIS 的裁剪表面反射率直方图，用于预冷泉火灾的所有波段。
从直方图中，您可以看到表面反射率的值范围似乎更合适，但仍不在 0 和 1 之间。
4、MODIS 影像中的反射率值
如博客前面所述，反射率值的正常范围是 0 到 1，其中 1 是最亮的值，0 是最暗的值。
再次查看上面为波段 1 计算的最小值和最大值。您注意到什么？
如您所见，最小值和最大值大大超出 0 到 1 的预期范围。查看波段 1 的直方图，您还可以看到值的范围不是您所期望的。
是什么原因造成的？ 要回答这个问题，您需要先更好地理解数据，然后才能更多地使用它。
（1）比例因子
使用遥感数据时，比例因子很常见。 数据很大，比例因子用于保持数据较小。 例如，存储带小数的数字（称为浮点数）比存储整数需要更多的空间。 因此，通常删除传感数据应用了一个比例因子，可用于
查看 MOD09GA 产品的表格，您可以看到 MODIS 数据的比例因子为 0.0001。 这意味着您应该将每一层乘以该值以获得数据的实际反射率值。
您可以使用 numpy 数组数学（有时在 GIS 工具中称为栅格数学）将此比例因子值应用于堆栈中的所有图层。 此处将整个数组乘以 .0001 以缩放每个图层或波段。
# Scale values of MODIS imagery stack
modis_bands_pre_scaled = modis_clip * 0.0001
# Identify minimum and maximum values of scaled band 1 (red)
print(modis_bands_pre_scaled[1].min(), modis_bands_pre_scaled[1].max())
<xarray.DataArray ()>
array(0.2496, dtype=float32)
Coordinates:
band int64 2
spatial_ref int64 0 <xarray.DataArray ()>
array(0.3013, dtype=float32)
Coordinates:
band int64 2
spatial_ref int64 0# Create a colors and titles list to use in the histogram
colors = ['r', 'k', 'b', 'g', 'k', 'y', 'y']
titles = ["Red Band", "Near Infrared (NIR) Band", "Blue/Green Band", "Green Band",
"Near Infrared (NIR) Band", "Mid-infrared Band", "Mid-infrared Band"]
# Plot histogram
ep.hist(modis_bands_pre_scaled.values,
colors=colors,
title=titles,
cols=2)
plt.show()
来自 MODIS 的裁剪和缩放表面反射率的直方图，用于预冷泉火灾的所有波段。
文章出处登录后可见！
