Process multiband rasters (Sentinel-with hndex and create indices-Tutoriel Python-php.cn

Hi, In previous blog we talked about how we can do raster analysis using h3 indexes and postgresql for a single band raster. In this blog we will talk about how we can process multiband raster and create indices with ease. We will be using sentinel-2 image and create NDVI from the processed h3 cells and visualize the results

Download sentinel 2 data

We are downloading the sentinel 2 data from https://apps.sentinel-hub.com/eo-browser/ in Pokhara, Nepal Area , Just to make sure lake is in the image grid so that it would be easy for us to validate the NDVI result

Process multiband rasters (Sentinel-with hndex and create indices

To download sentinel image with all bands :

You need to create an account
Find the image in your area select the grid that covers your Area of interest
Zoom to the grid, And click onicon on right vertical bar
After that go to analytical tab and select all the bands with image format as tiff 32 bit , high resolution , wgs1984 format and all bands checkd

Process multiband rasters (Sentinel-with hndex and create indices

You can also download pregenerated indices such as NDVI , False color tiff only or specific bands whichever best suits your need . We are downloading all the bands as we wanna do the processing by ourselves

Click download

Process multiband rasters (Sentinel-with hndex and create indices

Preprocess

We get all the bands as separate tiff from the sentinel as we downloaded raw format

Process multiband rasters (Sentinel-with hndex and create indices

lets create a composite image :

This can be done through GIS tools or gdal

Using gdal_merge:

We need to rename the downloaded file to band1,band2 like this to avoid slashes in the filename
Lets process upto band 9 for this exercise , you can choose the band as your requirement

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

Copier après la connexion

Using QGIS:

Load all individual bands to QGIS
Go to Raster > Miscellanaeous > Merge

Process multiband rasters (Sentinel-with hndex and create indices

While merging you need to make sure you check 'place each input file in sep band'

Process multiband rasters (Sentinel-with hndex and create indices

Now export your merged tiff to raw geotiff as composite

Housekeeping

Make sure your image is in WGS1984 in our case image is already is in ws1984 so no need for the conversion
Make sure you don't have any nodata if yes fill them with 0

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

Copier après la connexion

Finally make sure your ouput image is in COG

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

Copier après la connexion

I am using the bash script provided in cog2h3 repo to automate those

sudo bash pre.sh sentinel2_composite.tif

Copier après la connexion

Process and creation of h3 cells

Now finally as we have done the preprocessing script , lets move forward to compute h3 cells for each bands in the composite cog image

Install cog2h3

pip install cog2h3

Copier après la connexion

Export you database credentials

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

Copier après la connexion

We are using resolution 10 for this sentinel image, however you will also see in the script itself which will print the optimal resolution for your raster that makes the h3 cell smaller than your smallest of pixel in raster.

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

Copier après la connexion

It took a minute for us to compute and store result in postgresql

Logs :

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

Copier après la connexion

Analyze

Since now we have our data in postgresql , Lets do some analysis

Verify we have all the bands we processed ( Remember we processed from band 1 to 9 )

select * from sentinel

Copier après la connexion

Process multiband rasters (Sentinel-with hndex and create indices

Compute ndvi for each cell

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

Copier après la connexion

Query Plan :

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 |

Copier après la connexion

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 )

Copier après la connexion

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;

Copier après la connexion

Create index on geometry

create index on ndvi_sentinel(geometry);

Copier après la connexion

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

Process multiband rasters (Sentinel-with hndex and create indices

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 :

Process multiband rasters (Sentinel-with hndex and create indices

Check Lake

Process multiband rasters (Sentinel-with hndex and create indices
As there were clouds around the lake hence nearby fields are covered by cloud which makes sense