嗨，在之前的博客中，我们讨论了如何使用 h3 索引和 postgresql 对单波段栅格进行栅格分析。在本博客中，我们将讨论如何处理多波段栅格并轻松创建索引。我们将使用 Sentinel-2 图像并从处理后的 h3 细胞创建 NDVI 并可视化结果

下载哨兵2数据

我们正在从尼泊尔博卡拉地区的 https://apps.sentinel-hub.com/eo-browser/ 下载哨兵 2 数据，只是为了确保湖泊位于图像网格中，以便于我们验证NDVI结果

下载所有波段的哨兵图像：

• 您需要创建一个帐户
• 查找您所在区域的图像，选择覆盖您感兴趣区域的网格
• 缩放至网格，然后点击右侧竖条上的图标
• 之后转到分析选项卡并选择图像格式为 tiff 32 位、高分辨率、wgs1984 格式的所有波段并选中所有波段

您还可以下载预生成的指数，例如 NDVI、仅假色 tiff 或最适合您需要的特定波段。我们正在下载所有频段，因为我们想自己进行处理

• 点击下载

预处理

当我们下载原始格式时，我们将所有频段作为与哨兵分开的 tiff

• 让我们创建一个合成图像：

这可以通过 GIS 工具或 gdal 来完成

1. 使用 gdal_merge:

我们需要将下载的文件重命名为 band1,band2 以避免文件名中出现斜杠
本练习中我们最多处理频段 9，您可以根据需要选择频段

gdal_merge.py -separate -o sentinel2_composite.tif band1.tif band2.tif band3.tif band4.tif band5.tif band6.tif band7.tif band8.tif band9.tif 

1. 使用QGIS :

• 加载所有单独的波段到QGIS
• 转到栅格 > 其他 > 合并

• 合并时，您需要确保选中“将每个输入文件放入 sep band”

• 现在将合并的 tiff 作为复合材料导出到原始 geotiff

家政

• 确保您的图像采用 WGS1984 在我们的例子中，图像已经是 ws1984，所以不需要转换
• 确保您没有任何 nodata 如果有则用 0 填充

  gdalwarp -overwrite -dstnodata 0 "$input_file" "${output_file}_nodata.tif"

• 最后确保您的输出图像是 COG 格式

  gdal_translate -of COG "$input_file" "$output_file"

我正在使用 cog2h3 repo 中提供的 bash 脚本来自动化这些

sudo bash pre.sh sentinel2_composite.tif

h3细胞的处理和创建

现在，我们终于完成了预处理脚本，让我们继续计算复合齿轮图像中每个波段的 h3 单元格

• 安装cog2h3

  pip install cog2h3

• 导出您的数据库凭据

  export DATABASE_URL="postgresql://user:password@host:port/database"

• 跑步

我们对此哨兵图像使用分辨率 10，但是您还会在脚本本身中看到，它将打印栅格的最佳分辨率，使 h3 单元小于栅格中的最小像素。

  cog2h3 --cog sentinel2_composite_preprocessed.tif --table sentinel --multiband --res 10

我们花了一分钟的时间来计算结果并将其存储在 postgresql 中

日志：

2024-08-24 08:39:43,233 - INFO - Starting processing
2024-08-24 08:39:43,234 - INFO - COG file already exists at sentinel2_composite_preprocessed.tif
2024-08-24 08:39:43,234 - INFO - Processing raster file: sentinel2_composite_preprocessed.tif
2024-08-24 08:39:43,864 - INFO - Determined Min fitting H3 resolution for band 1: 11
2024-08-24 08:39:43,865 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:44,037 - INFO - Resampling Done for band 1
2024-08-24 08:39:44,037 - INFO - New Native H3 resolution for band 1: 10
2024-08-24 08:39:44,738 - INFO - Calculation done for res:10 band:1
2024-08-24 08:39:44,749 - INFO - Determined Min fitting H3 resolution for band 2: 11
2024-08-24 08:39:44,749 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:44,757 - INFO - Resampling Done for band 2
2024-08-24 08:39:44,757 - INFO - New Native H3 resolution for band 2: 10
2024-08-24 08:39:45,359 - INFO - Calculation done for res:10 band:2
2024-08-24 08:39:45,366 - INFO - Determined Min fitting H3 resolution for band 3: 11
2024-08-24 08:39:45,366 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:45,374 - INFO - Resampling Done for band 3
2024-08-24 08:39:45,374 - INFO - New Native H3 resolution for band 3: 10
2024-08-24 08:39:45,986 - INFO - Calculation done for res:10 band:3
2024-08-24 08:39:45,994 - INFO - Determined Min fitting H3 resolution for band 4: 11
2024-08-24 08:39:45,994 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:46,003 - INFO - Resampling Done for band 4
2024-08-24 08:39:46,003 - INFO - New Native H3 resolution for band 4: 10
2024-08-24 08:39:46,605 - INFO - Calculation done for res:10 band:4
2024-08-24 08:39:46,612 - INFO - Determined Min fitting H3 resolution for band 5: 11
2024-08-24 08:39:46,612 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:46,619 - INFO - Resampling Done for band 5
2024-08-24 08:39:46,619 - INFO - New Native H3 resolution for band 5: 10
2024-08-24 08:39:47,223 - INFO - Calculation done for res:10 band:5
2024-08-24 08:39:47,230 - INFO - Determined Min fitting H3 resolution for band 6: 11
2024-08-24 08:39:47,230 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:47,239 - INFO - Resampling Done for band 6
2024-08-24 08:39:47,239 - INFO - New Native H3 resolution for band 6: 10
2024-08-24 08:39:47,829 - INFO - Calculation done for res:10 band:6
2024-08-24 08:39:47,837 - INFO - Determined Min fitting H3 resolution for band 7: 11
2024-08-24 08:39:47,837 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:47,845 - INFO - Resampling Done for band 7
2024-08-24 08:39:47,845 - INFO - New Native H3 resolution for band 7: 10
2024-08-24 08:39:48,445 - INFO - Calculation done for res:10 band:7
2024-08-24 08:39:48,453 - INFO - Determined Min fitting H3 resolution for band 8: 11
2024-08-24 08:39:48,453 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:48,461 - INFO - Resampling Done for band 8
2024-08-24 08:39:48,461 - INFO - New Native H3 resolution for band 8: 10
2024-08-24 08:39:49,046 - INFO - Calculation done for res:10 band:8
2024-08-24 08:39:49,054 - INFO - Determined Min fitting H3 resolution for band 9: 11
2024-08-24 08:39:49,054 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:49,062 - INFO - Resampling Done for band 9
2024-08-24 08:39:49,063 - INFO - New Native H3 resolution for band 9: 10
2024-08-24 08:39:49,647 - INFO - Calculation done for res:10 band:9
2024-08-24 08:39:51,435 - INFO - Converting H3 indices to hex strings
2024-08-24 08:39:51,906 - INFO - Overall raster calculation done in 8 seconds
2024-08-24 08:39:51,906 - INFO - Creating or replacing table sentinel in database
2024-08-24 08:40:03,153 - INFO - Table sentinel created or updated successfully in 11.25 seconds.
2024-08-24 08:40:03,360 - INFO - Processing completed

分析

现在我们的数据已经在postgresql中了，让我们做一些分析

• 验证我们是否拥有已处理的所有频段（请记住我们从频段 1 到 9 进行了处理）

select *
from sentinel

• 计算每个单元的 ndvi

explain analyze 
select h3_ix , (band8-band4)/(band8 band4) as ndvi
from public.sentinel

查询计划：

QUERY PLAN                                                                                                       |
----------------------------------------------------------------------------------------------------------------- 
Seq Scan on sentinel  (cost=0.00..28475.41 rows=923509 width=16) (actual time=0.014..155.049 rows=923509 loops=1)|
Planning Time: 0.080 ms                                                                                          |
Execution Time: 183.764 ms                                                                                       |

正如您在此处看到的那样，对于该区域中的所有行，计算都是即时的。对于所有其他索引都是如此，您可以使用 h3_ix 主键计算与其他表的复杂索引连接，并从中导出有意义的结果，而不必担心，因为 postgresql 能够处理复杂的查询和表连接。

可视化和验证

让我们可视化并验证计算的索引是否正确

• 创建表格（用于在 QGIS 中可视化）

create table ndvi_sentinel
as(
select h3_ix , (band8-band4)/(band8 band4) as ndvi
from public.sentinel )

• 让我们添加几何图形来可视化 h3 单元格 这仅是在 QGIS 中可视化所必需的，如果您自己构建一个最小的 API，则不需要它，因为您可以直接从查询构造几何图形

ALTER TABLE ndvi_sentinel  
ADD COLUMN geometry geometry(Polygon, 4326) 
GENERATED ALWAYS AS (h3_cell_to_boundary_geometry(h3_ix)) STORED;

• 在几何体上创建索引

create index on ndvi_sentinel(geometry);

• 在QGIS中连接数据库并根据ndvi值可视化表格 让我们获取费瓦湖或云附近的区域

据我们所知，-1.0 到 0.1 之间的值应该代表深水或浓密的云层
让我们看看这是否属实（使第一个类别透明以查看底层图像）

• 检查云：

• 检查湖

由于湖周围有云，因此附近的田野被云覆盖，这是有道理的

感谢您的阅读！下一篇博客见