摘要: 可视化栅格数据对于评估分析和模型的结果非常重要。可视化栅格数据的最佳方法之一是将其绘制在地图上，可在地图上与栅格数据和其他相关图层进行交互。 在 Python 中执行分析时，创建交互式地图可能很困难。这通常意味着必须从 Python 脚本或程序导出数据图层...
可视化栅格数据对于评估分析和模型的结果非常重要。可视化栅格数据的最佳方法之一是将其绘制在地图上，可在地图上与栅格数据和其他相关图层进行交互。
在 Python 中执行分析时，创建交互式地图可能很困难。这通常意味着必须从 Python 脚本或程序导出数据图层并在其他 GIS 软件（例如 ArcGIS 或 QGIS）中打开才能查看。
Leafmap是一个 Python 包，可直接在 Jupyter Notebook 中创建交互式地图。它允许添加保存到磁盘和 xarray 数据的栅格图层。Leafmap甚至支持分屏滑动视图以轻松比较不同的数据层。
接下来可按照以下步骤使用 leafmap 创建具有分屏、滑动视图的交互式 Python 地图。
数据
在以下内容中，将使用 Sentinel 2A 中的图像，以真彩色和假彩色显示该图像，同时还将使用 Sentinel 2A 数据来计算 NDVI。
接下来需要将 Sentinel 2A 数据下载到您的计算机上才能继续操作。哨兵图像可以从各种来源免费下载，不过很多人更喜欢直接访问源头，哥白尼浏览器（详情参见：https://browser.dataspace.copernicus.eu/）。
环境和包
需要一些包来执行本教程中的可视化和分析。建议使用minconda 发行版（详情参见：https://opensourceoptions.com/the-easiest-way-to-install-python-on-windows-11-also-works-for-mac-and-linux/），并使用 conda 和 mamba 安装软件包。
下面的代码块列出了要安装和导入的包，
from osgeo import gdal
import leafmap
import xarray as xr
import rioxarray
C:\Users\Konrad\miniconda3\envs\split\lib\site-packages\pandas\core\computation\expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.7.3' currently installed).
from pandas.core.computation.check import NUMEXPR_INSTALLED
创建虚拟栅格
Sentinel 收集 13 个不同频段的数据。每个波段包含不同波长光的反射率数据，这些数据存储在单独的光栅文件中。
为了显示哨兵图像，我们希望将它们包含在同一个文件中。因此将使用虚拟光栅格式来完成此操作。虚拟格式基本上是一个元数据文件，它从多个栅格获取信息，以便可以将它们组合起来，而无需实际处理底层数据。
我们将仅使用 2-5 波段（蓝、绿、红、近红外）创建一个虚拟栅格。
# Create a virtual raster from Sentinel 2a image
# path to a Sentinel 2A image
sen2a_base = r"path\to\sentinel\image\S2A_MSIL1C_20230626T180921_N0509_R084_T12TVK_20230626T215620.SAFE\S2A_MSIL1C_20230626T180921_N0509_R084_T12TVK_20230626T215620.SAFE\GRANULE\L1C_T12TVK_A041837_20230626T181453\IMG_DATA\T12TVK_20230626T180921_B{0}.jp2"
sen2a_fns = []
vrt_fn = "../data/my.vrt"
for band_num in range(2, 6):
sen2a_fns.append(sen2a_base.format(str(band_num).zfill(2)))
sen2_vrt = gdal.BuildVRT(vrt_fn, sen2a_fns, separate=True)
sen2_vrt = None
一旦 VRT 构建完成，将变量设置为None文件写入磁盘。
计算植被指数
接下来计算 NDVI 作为xarray数据数组。我们可以显示这种数据类型，而不需要将文件写入磁盘。
使用 rioxarray 将VRT作为 xarray 打开。
应用波段 3（红色）和 4（近红外）的 NDVI 方程，如下所示。
sen2_xr = rioxarray 。open_rasterio ( vrt_fn )
ndvi = ( sen2_xr . sel ( band = 4 ) - sen2_xr . sel ( band = 3 )) / ( sen2_xr . sel ( band = 4 ) + sen2_xr . sel ( band = 3 ))
现在使用 leafmap 创建地图。
首先创建一个 leafmap.Map 对象，然后应用该 split_map 函数。现在指定左层 ( ndvi )、右层（输入 VRT 的文件名）、左标签和右标签。左右参数是与 leafmap.add_raster_layer 相对应的参数字典。
对于哨兵图像，指定要应用的标注栏组合。[3，2，1]的组合将给出真实的颜色，[4，3，2]将给出彩色红外。运行代码，您的地图就会准备好, 不过加载地图和数据可能需要几分钟时间。
m = leafmap.Map()
m.split_map(
ndvi,
vrt_fn,
left_label="NDVI",
right_label="Sentinel 2A",
left_args={
"layer_name": "NDVI",
"colormap": "turbo"
},
right_args={
"layer_name": "Sentinel 2A",
"bands": [4, 3, 2]
}
)
m
Map(center=[20, 0], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title', 'zoom_out_text…
总结
Leafmap 为地理空间数据可视化提供了一个优秀的 Python 包。与配合使用时xarrary，它还提供了一种查看分析结果的好方法，而无需将数据写入磁盘并在不同的应用程序中打开它。
通过滑动的分割地图为数据可视化和比较提供了良好的触感。
