嗨，在之前的博客中，我们讨论了如何使用 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来完成
使用 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
使用 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 processing2024-08-24 08:39:43,234 - info - cog file already exists at sentinel2_composite_preprocessed.tif2024-08-24 08:39:43,234 - info - processing raster file: sentinel2_composite_preprocessed.tif2024-08-24 08:39:43,864 - info - determined min fitting h3 resolution for band 1: 112024-08-24 08:39:43,865 - info - resampling original raster to: 200.786148m2024-08-24 08:39:44,037 - info - resampling done for band 12024-08-24 08:39:44,037 - info - new native h3 resolution for band 1: 102024-08-24 08:39:44,738 - info - calculation done for res:10 band:12024-08-24 08:39:44,749 - info - determined min fitting h3 resolution for band 2: 112024-08-24 08:39:44,749 - info - resampling original raster to: 200.786148m2024-08-24 08:39:44,757 - info - resampling done for band 22024-08-24 08:39:44,757 - info - new native h3 resolution for band 2: 102024-08-24 08:39:45,359 - info - calculation done for res:10 band:22024-08-24 08:39:45,366 - info - determined min fitting h3 resolution for band 3: 112024-08-24 08:39:45,366 - info - resampling original raster to: 200.786148m2024-08-24 08:39:45,374 - info - resampling done for band 32024-08-24 08:39:45,374 - info - new native h3 resolution for band 3: 102024-08-24 08:39:45,986 - info - calculation done for res:10 band:32024-08-24 08:39:45,994 - info - determined min fitting h3 resolution for band 4: 112024-08-24 08:39:45,994 - info - resampling original raster to: 200.786148m2024-08-24 08:39:46,003 - info - resampling done for band 42024-08-24 08:39:46,003 - info - new native h3 resolution for band 4: 102024-08-24 08:39:46,605 - info - calculation done for res:10 band:42024-08-24 08:39:46,612 - info - determined min fitting h3 resolution for band 5: 112024-08-24 08:39:46,612 - info - resampling original raster to: 200.786148m2024-08-24 08:39:46,619 - info - resampling done for band 52024-08-24 08:39:46,619 - info - new native h3 resolution for band 5: 102024-08-24 08:39:47,223 - info - calculation done for res:10 band:52024-08-24 08:39:47,230 - info - determined min fitting h3 resolution for band 6: 112024-08-24 08:39:47,230 - info - resampling original raster to: 200.786148m2024-08-24 08:39:47,239 - info - resampling done for band 62024-08-24 08:39:47,239 - info - new native h3 resolution for band 6: 102024-08-24 08:39:47,829 - info - calculation done for res:10 band:62024-08-24 08:39:47,837 - info - determined min fitting h3 resolution for band 7: 112024-08-24 08:39:47,837 - info - resampling original raster to: 200.786148m2024-08-24 08:39:47,845 - info - resampling done for band 72024-08-24 08:39:47,845 - info - new native h3 resolution for band 7: 102024-08-24 08:39:48,445 - info - calculation done for res:10 band:72024-08-24 08:39:48,453 - info - determined min fitting h3 resolution for band 8: 112024-08-24 08:39:48,453 - info - resampling original raster to: 200.786148m2024-08-24 08:39:48,461 - info - resampling done for band 82024-08-24 08:39:48,461 - info - new native h3 resolution for band 8: 102024-08-24 08:39:49,046 - info - calculation done for res:10 band:82024-08-24 08:39:49,054 - info - determined min fitting h3 resolution for band 9: 112024-08-24 08:39:49,054 - info - resampling original raster to: 200.786148m2024-08-24 08:39:49,062 - info - resampling done for band 92024-08-24 08:39:49,063 - info - new native h3 resolution for band 9: 102024-08-24 08:39:49,647 - info - calculation done for res:10 band:92024-08-24 08:39:51,435 - info - converting h3 indices to hex strings2024-08-24 08:39:51,906 - info - overall raster calculation done in 8 seconds2024-08-24 08:39:51,906 - info - creating or replacing table sentinel in database2024-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 ndvifrom 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_sentinelas(select h3_ix , (band8-band4)/(band8+band4) as ndvifrom 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 之间的值应该代表深水或浓密的云层
让我们看看这是否属实（使第一个类别透明以查看底层图像）
检查云：
检查湖
由于湖周围有云，因此附近的田野被云覆盖，这是有道理的
感谢您的阅读！下一篇博客见
以上就是处理多波段栅格（Sentinel-使用 hndex 并创建索引的详细内容，更多请关注创想鸟其它相关文章！
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。
如发现本站有涉嫌抄袭侵权/违法违规的内容， 请发送邮件至 chuangxiangniao@163.com 举报，一经查实，本站将立刻删除。
发布者：程序猿，转转请注明出处：https://www.chuangxiangniao.com/p/1349465.html
微信扫一扫 
支付宝扫一扫 
