Verarbeiten Sie Multiband-Raster (Sentinel-mit hndex und erstellen Sie Indizes

WBOY
Freigeben: 2024-08-25 06:01:33
Original
900 Leute haben es durchsucht

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 on Process multiband rasters (Sentinel-with hndex and create indices icon 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

  1. 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 
Nach dem Login kopieren
  1. 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"
Nach dem Login kopieren
  • Finally make sure your ouput image is in COG
  gdal_translate -of COG "$input_file" "$output_file"
Nach dem Login kopieren

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

sudo bash pre.sh sentinel2_composite.tif
Nach dem Login kopieren

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
Nach dem Login kopieren
  • Export you database credentials
  export DATABASE_URL="postgresql://user:password@host:port/database"
Nach dem Login kopieren
  • Run

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
Nach dem Login kopieren

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
Nach dem Login kopieren

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
Nach dem Login kopieren

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
Nach dem Login kopieren

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                                                                                       |
Nach dem Login kopieren

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 )
Nach dem Login kopieren
  • 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;
Nach dem Login kopieren
  • Create index on geometry
create index on ndvi_sentinel(geometry);
Nach dem Login kopieren
  • 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

Process multiband rasters (Sentinel-with hndex and create indices

Thank you for reading ! See you in next blog

Das obige ist der detaillierte Inhalt vonVerarbeiten Sie Multiband-Raster (Sentinel-mit hndex und erstellen Sie Indizes. Für weitere Informationen folgen Sie bitte anderen verwandten Artikeln auf der PHP chinesischen Website!

Quelle:dev.to
Erklärung dieser Website
Der Inhalt dieses Artikels wird freiwillig von Internetnutzern beigesteuert und das Urheberrecht liegt beim ursprünglichen Autor. Diese Website übernimmt keine entsprechende rechtliche Verantwortung. Wenn Sie Inhalte finden, bei denen der Verdacht eines Plagiats oder einer Rechtsverletzung besteht, wenden Sie sich bitte an admin@php.cn
Beliebte Tutorials
Mehr>
Neueste Downloads
Mehr>
Web-Effekte
Quellcode der Website
Website-Materialien
Frontend-Vorlage