Привет! В предыдущем блоге мы говорили о том, как можно выполнить растровый анализ с использованием индексов h3 и postgresql для одноканального растра. В этом блоге мы поговорим о том, как можно легко обрабатывать многоканальные растры и создавать индексы. Мы будем использовать изображение Sentinel-2, создадим NDVI из обработанных ячеек h3 и визуализируем результаты

Загрузить данные Sentinel 2

Мы загружаем данные Sentinel 2 с https://apps.sentinel-hub.com/eo-browser/ в Покхаре, Непал. Просто чтобы убедиться, что озеро находится в сетке изображений, чтобы его можно было легко нам для проверки результата NDVI

Чтобы загрузить изображение дозорного со всеми диапазонами:

• Вам необходимо создать учетную запись
• Найдите изображение в своей области, выберите сетку, охватывающую вашу область интересов
• Приблизьтесь к сетке и нажмите значок на правой вертикальной панели
• После этого перейдите на вкладку «Аналитика» и выберите все каналы с форматом изображения tiff 32 бит, высоким разрешением, форматом wgs1984 и все каналы отмечены галочкой

Вы также можете загрузить предварительно созданные индексы, такие как NDVI, только Tiff цвета False или определенные полосы, в зависимости от того, что лучше всего соответствует вашим потребностям. Мы скачиваем все группы, так как хотим сами выполнить обработку

• Нажмите «Загрузить»

Предварительная обработка

Мы получаем все группы в виде отдельного TIFF от Sentinel, поскольку мы загружали необработанный формат

• давайте создадим составное изображение:

Это можно сделать с помощью инструментов ГИС или 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
• Перейдите в раздел «Растр» > «Разное» > «Объединить»

• При объединении необходимо убедиться, что вы установили флажок «помещать каждый входной файл в полосу разделения»

• Теперь экспортируйте объединенный TIFF в необработанный геотифф как составной

Уборка номеров

• Убедитесь, что ваше изображение находится в формате WGS1984. в нашем случае изображение уже есть в ws1984, поэтому преобразование не требуется
• Убедитесь, что у вас нет никаких данных, если да, введите 0

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

• Наконец, убедитесь, что ваше выходное изображение находится в COG.

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

Для автоматизации этих операций я использую скрипт bash из репозитория cog2h3.

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 должно обозначать глубокую воду или плотные облака.
давайте посмотрим, так ли это (сделаем первую категорию прозрачной, чтобы увидеть основное изображение)

• Проверить облака:

• Чек-Лейк

Поскольку вокруг озера были облака, близлежащие поля были покрыты облаками, что имеет смысл

Спасибо, что читаете! Увидимся в следующем блоге