處理多波段柵格（Sentinel-with hndex 並建立索引
發佈於2024-08-25
瀏覽：498
嗨，在之前的博客中，我们讨论了如何使用 h3 索引和 postgresql 对单波段栅格进行栅格分析。在本博客中，我们将讨论如何处理多波段栅格并轻松创建索引。我们将使用 Sentinel-2 图像并从处理后的 h3 细胞创建 NDVI 并可视化结果
下载哨兵2数据
我们正在从尼泊尔博卡拉地区的 https://apps.sentinel-hub.com/eo-browser/ 下载哨兵 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 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 之间的值应该代表深水或浓密的云层
让我们看看这是否属实（使第一个类别透明以查看底层图像）
- 检查云：
- 检查湖
由于湖周围有云，因此附近的田野被云覆盖，这是有道理的
感谢您的阅读！下一篇博客见
版本聲明
本文轉載於：https://dev.to/krschap/process-multiband-rasters-sentinel-2-with-h3-index-and-create-indices-bem?1如有侵犯，請聯絡[email protected]刪除
最新教學
更多>
-
如何使用Regex在PHP中有效地提取括號內的文本php：在括號內提取文本在處理括號內的文本時，找到最有效的解決方案是必不可少的。一種方法是利用PHP的字符串操作函數，如下所示： 作為替代 $ text ='忽略除此之外的一切（text）'; preg_match（'＃（（。 &&& [Regex使用模式來搜索特...程式設計 發佈於2025-04-12 -
如何將MySQL數據庫添加到Visual Studio 2012中的數據源對話框中？在Visual Studio 2012 儘管已安裝了MySQL Connector v.6.5.4，但無法將MySQL數據庫添加到實體框架的“ DataSource對話框”中。為了解決這一問題，至關重要的是要了解MySQL連接器v.6.5.5及以後的6.6.x版本將提供MySQL的官方Visual...程式設計 發佈於2025-04-12
-
如何檢查對像是否具有Python中的特定屬性？方法來確定對象屬性存在尋求一種方法來驗證對像中特定屬性的存在。考慮以下示例，其中嘗試訪問不確定屬性會引起錯誤： >>> a = someClass（） >>> A.property Trackback（最近的最新電話）： 文件“ ”，第1行， AttributeError: SomeClass...程式設計 發佈於2025-04-12
-
使用CSS3和SVG創建帶邊框的波浪形狀在嘗試使用形狀使用形狀的CSS3設計CSS3時，在CSS3 中創建一個波形形狀，由於所需結果可能無法獲得邊框和背景彩色設置的限制，因此無法實現所需的結果。要克服這一點，請考慮使用SVG代替DIV來進行波形。 實現：在容器中，將內容和SVG放置在波形中。將SVG向右浮動。 svg設計： 很大，並用...程式設計 發佈於2025-04-12
-
Android如何向PHP服務器發送POST數據？在android apache httpclient（已棄用） httpclient httpclient = new defaulthttpclient（）; httppost httppost = new httppost（“ http://www.yoursite.com/script.p...程式設計 發佈於2025-04-12
-
版本5.6.5之前，使用current_timestamp與時間戳列的current_timestamp與時間戳列有什麼限制？在時間戳列上使用current_timestamp或MySQL版本中的current_timestamp或在5.6.5 此限制源於遺留實現的關注，這些限制需要對當前的_timestamp功能進行特定的實現。 創建表`foo`（ `Productid` int（10）unsigned not ...程式設計 發佈於2025-04-12
-
如何將多種用戶類型（學生，老師和管理員）重定向到Firebase應用中的各自活動？Red: How to Redirect Multiple User Types to Respective ActivitiesUnderstanding the ProblemIn a Firebase-based voting app with three distinct user type...程式設計 發佈於2025-04-12
-
MySQL中基數如何影響索引優化？在mySQL cardinital可以分為兩個類別：高和低。具有較高基數的列具有大量唯一值，而低心電圖列的不同值數量有限。 Cardinality and Index OptimizationCardinality is closely related to indexing, whic...程式設計 發佈於2025-04-12
-
前端挑戰：提陞技能的實戰指南提升前端開發技能的最佳途徑？那就是實踐！動手構建網站，這是最有效的學習方法。如果能從中獲得報酬，那就再好不過了；即使是為自身或親友構建網站，也能顯著提陞技能。即使只是為了練習而創建項目，也能讓你快速成長，這絕對比單純閱讀資料有效得多！ 以下是一些資源，它們鼓勵你通過構建項目來提陞技能： Fron...程式設計 發佈於2025-04-12
-
PostgreSQL中如何提取每個ID的最後一行數據？在postgresql To accomplish this in Postgresql, two methods are commonly used:Distinct On OperatorPostgresql provides the distinct on operator, which ...程式設計 發佈於2025-04-12
-
如何在其容器中為DIV創建平滑的左右CSS動畫？通用CSS動畫，用於左右運動 ，我們將探索創建一個通用的CSS動畫，以向左和右移動DIV，從而到達其容器的邊緣。該動畫可以應用於具有絕對定位的任何div，無論其未知長度如何。 問題：使用左直接導致瞬時消失 更加流暢的解決方案：混合轉換和左 [並實現平穩的，線性的運動，我們介紹了線性的轉換。...程式設計 發佈於2025-04-12
-
您可以使用CSS在Chrome和Firefox中染色控制台輸出嗎？在javascript console 中顯示顏色是可以使用chrome的控制台顯示彩色文本，例如紅色的redors，for for for for錯誤消息？ 回答是的，可以使用CSS將顏色添加到Chrome和Firefox中的控制台顯示的消息（版本31或更高版本）中。要實現這一目標，請使用以下...程式設計 發佈於2025-04-12
-
如何使用PHP將斑點（圖像）正確插入MySQL？essue VALUES('$this->image_id','file_get_contents($tmp_image)')";This code builds a string in PHP, but the function call fil...程式設計 發佈於2025-04-12
-
IE6兼容自定義數據屬性指南在IE 6中的自定義數據屬性：分發神話 custom Data屬性，HTML5的關鍵功能，使開發人員能夠將非可視數據附加到HTML elements for HTML elements以供以後回顧。但是，一個常見的誤解包圍著它們與諸如Internet Explorer 6的傳統瀏覽器的兼容性。 ...程式設計 發佈於2025-04-12
-
如何同步迭代並從PHP中的兩個等級陣列打印值？同步的迭代和打印值來自相同大小的兩個數組使用兩個數組相等大小的selectbox時，一個包含country代碼的數組，另一個包含鄉村代碼，另一個包含其相應名稱的數組，可能會因不當提供了exply for for for the uncore for the forsion for for ytry...程式設計 發佈於2025-04-12
學習中文
- 1 走路用中文怎麼說？ 走路中文發音，走路中文學習
- 2 坐飛機用中文怎麼說？ 坐飞机中文發音，坐飞机中文學習
- 3 坐火車用中文怎麼說？ 坐火车中文發音，坐火车中文學習
- 4 坐車用中文怎麼說？ 坐车中文發音，坐车中文學習
- 5 開車用中文怎麼說？ 开车中文發音，开车中文學習
- 6 游泳用中文怎麼說？ 游泳中文發音，游泳中文學習
- 7 騎自行車用中文怎麼說？ 骑自行车中文發音，骑自行车中文學習
- 8 你好用中文怎麼說？ 你好中文發音，你好中文學習
- 9 謝謝用中文怎麼說？ 谢谢中文發音，谢谢中文學習
- 10 How to say goodbye in Chinese? 再见Chinese pronunciation, 再见Chinese learning
