处理多波段栅格（Sentinel-使用 hndex 并创建索引

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

下载哨兵2数据

我们正在从尼泊尔博卡拉地区的https://apps.sentinel-hub.com/eo-browser/下载sentinel 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 存储库中提供的 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                                                                                       |

As you can see here for all the rows in that area the calculation is instant . This is true for all other indices and you can compute complex indices join with other tables using the h3_ix primary key and derive meaningful result out of it without worrying as postgresql is capable of handling complex queries and table join.

Visualize and verification

Lets visualize and verify if the computed indices are true

• Create table ( for visualizing in QGIS )

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

• Lets add geometry to visualize the h3 cells This is only necessary to visualize in QGIS , if you build an minimal API by yourself you don't need this as you can construct geometry directly from query

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

• Create index on geometry

create index on ndvi_sentinel(geometry);

• Connect your database in QGIS and visualize the table on the basis of ndvi value Lets get the area near Fewa lake or cloud

As we know value between -1.0 to 0.1 should represent Deep water or dense clouds
lets see if thats true ( making first category as transparent to see the underlying image )

• Check clouds :

• Check Lake

As there were clouds around the lake hence nearby fields are covered by cloud which makes sense

Thank you for reading ! See you in next blog

以上是处理多波段栅格（Sentinel-使用 hndex 并创建索引的详细内容。更多信息请关注PHP中文网其他相关文章！