可视化栅格数据对于评估分析和模型的结果非常重要。 可视化栅格数据的最佳方法之一是将其绘制在地图上， 可在地图上与栅格数据和其他相关图层进行交互。
在 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 安装软件包。
下面的代码块列出了要安装和导入的包，
<pre> from osgeo import gdal import leafmap import xarray as xr import rioxarray </pre>
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]将给出彩色红外。 运行代码，您的地图就会准备好, 不过加载地图和数据可能需要几分钟时间。
<pre> 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 </pre>
Map(center=[20, 0], controls=(ZoomControl(options=['position', 'zoom_in_text', 'zoom_in_title', 'zoom_out_text…
总结
Leafmap 为地理空间数据可视化提供了一个优秀的 Python 包。 与配合使用时xarrary， 它还提供了一种查看分析结果的好方法， 而无需将数据写入磁盘并在不同的应用程序中打开它。
通过滑动的分割地图为数据可视化和比较提供了良好的触感。