SWOT 水资源卫星点云数据处理

表层水和海洋地形学任务 Surface Water Ocean Topography mission (SWOT) ，下文称之为SWOT

该任务于 2022 年 12 月 16 日启动，有望改变水位（海洋和内陆）监测的突破性技术。该任务首次通过提供比以往任何时候都更精确的水位高度测量来研究地球表面的几乎所有水。

该任务于 2022 年推出，首次公开版本已于 2023 年 12 月 5 日宣布（详情见：https://swot.jpl.nasa.gov/news/113/swot-data-first-public-release/），距今马上一年，虽然处于“非常早期的阶段”并且具有“已知限制”（如发行说明）

马上我们将有一年满的数据，到目前为止，还没有看到对这些首批发布的产品的评估论文。

由于种种原因，国内水文学者对SWOT数据认识也不足，因此，本文章提供这些数据的快速教程。

什么是SWOT卫星？

快问快答

• SWOT卫星如何运作的？

该任务的目标是利用两种不同形式的卫星技术（即雷达测高和激光雷达图像）创建世界表层水的地图。

• 都有什么数据？

在雷达测高法中，使用雷达测量高程随时间的变化，而在激光雷达成像中，使用激光测量地球表面物体之间的距离。这些技术结合起来，将产生我们星球上海洋、河流、湖泊和湿地的高精度图片。

• 有什么潜在应用？

以更好地了解不同水体如何相互作用，以及它们所处的环境如何。也可以使用它在洪水或干旱发生之前进行预测。也可能有助于缓解一些地区因气候变化（如厄尔尼诺事件或极地冰川融化）的影响而造成的饮用水短缺或干旱状况。

这项任务有可能大大提高与现在和未来负责任地管理全球水资源相关的科学理解和实际应用，因为它将让我们从太空全面了解地球的水文循环。

SWOT数据简介

数据以点云的形式储存在.nc文件中

数据包括：高频数据的初级产品（L1B_HR_SLC）和基于L1B_HR_SLC生产的二级产品（L2_HR_）

二级产品仅保留了陆地内的水体部分，因此可以用来研究精细尺度的水文现象。「对于绝大多数水文学者，二级产品就足够了。」

二级产品中提供了河流水面高度、宽度、坡度、流量等信息；也有湖泊水位、面积、水深、水储量等信息

有两种主要的数据格式，一种是点云，如L2_HR_PIXC；另一种是河流，如L2_HR_RiverSP

这次先介绍如何处理点云数据，下次我们再详细介绍如何处理河流数据

数据读取与下载

由于数据量很大，所有数据本地下载也不太现实。第一步是定义感兴趣区域 （AOI） 进行分析（这里用鄱阳湖为例）。

https://geojson.io/

推荐用geojson选择感兴趣的区域

然后我们用earthaccess 包直接下载数据处理，这需要先注册。

from pathlib import Path

# Enter username and password to login
earthaccess.login()

# Search for products within the AOI and Time frame
results = earthaccess.search_data(short_name='SWOT_L2_HR_PIXC_2.0', 
                                  temporal=('2023-08-01', '2023-08-30'), 
                                  bounding_box=aoi.bounds)

# Display the granules found
items = [item['meta']['native-id'] for item in results]
earthaccess.download(results[:1], 'D:/Onedrive/Untitled Folder/python/SWOT/data/swot')

这样数据就下载好了

数据量还挺大的，仅仅一个湖就快1G了

接下来我们提取所需的字段，如经纬度，高度，水体分类等等，将其转为geodataframe

部分代码如下：

ds_pixc = xr.open_dataset('SWOT_L2_HR_PIXC_001_534_103L_20230809T065206_20230809T065217_PGC0_01.nc')

# Load a GeoPandas DataFrame with the points and some variables
gdf = gpd.GeoDataFrame(
    data={
        'height': ds_pixc.height.values,
        'classification': ds_pixc.classification.values
        xxxx加上其它你需要的字段等等
    },
    )
)

数据绘制

我们想看看点云数据是如何的，代码如下：

custom_colors = ListedColormap([ 'navy', 'green', 'darkgreen', 'red', 'orange', 'blue', 'royalblue'])

classes = {
    1: 'land',
    2: 'land_near_water',
    3: 'water_near_land',
    4: 'open_water',
    5: 'dark_water'
}

gdf['classes'] = gdf['classification'].map(classes).astype('category')

# apply Log to the coherent_power to better visualization
gdf['log_power'] = np.log(gdf['coherent_power'])

gdf.plot(markersize=0.1, column='classes', figsize=(10, 10), legend=True, cmap=custom_colors)

可以看到数据的精度还是挺高的，但其实点云数据是不太好分析的，我们还是得把他转成tif才好操作

数据转换

首先获取裁剪后点云数据的边界，然后设定栅格分辨率（这里设置为30米），为栅格化创建形状和对应值的生成器，最后生成TIFF

由于篇幅原因，只展示核心代码：

with rasterio.open(
    'sampled_point_cloud_10m.tif',
    'w',
    driver='GTiff',
    height=height,
    width=width,
    count=1,
    dtype=rasterio.uint8,
    crs="EPSG:4326",
    transform=transform,
) as dst:
    # 使用rasterize函数将点云分类写入栅格
    raster = rasterize(
        shapes=shapes,
        out_shape=(height, width),
        transform=transform,
        fill=0,  # 设置背景值
        dtype=rasterio.uint8
    )
    
    # 写入栅格数据
    dst.write(raster, 1)

print("10米分辨率的GeoTIFF已成功生成并保存为 'sampled_point_cloud_30m.tif'")

最后输出了栅格，我们打开QGIS看一下效果怎么样

这样我们就把他成功栅格化了，重复上述操作，我们就能生成一张全球的高精度水体分布图，从而正确地分析和研究水资源变化。