Skip to content
GitLab
Commit 2ca9138c authored by Ose Kenji's avatar Ose Kenji
Browse files

use of dinamis-sdk, from image search to ndvi

parent a863a649
main
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1250 additions and 0 deletions
+1250 -0
examples/dinamis_sdk_ndvi.ipynb 0 → 100644
+ 1250
− 0
View file @ 2ca9138c
%% Cell type:markdown id:99c09edf-cdb9-4044-8579-736a14138e0c tags:
# Exemple d'utilisation PySTAC Client avec DINAMIS-SDK
---
auteur : Kenji Ose, UMR-TETIS, INRAE
date : 2023-06-16
---
%% Cell type:markdown id:b1863dc0-fe72-44f5-bc54-b5e4c8d3b32c tags:
## 1. Dinamis-SDK
### 1.1. Brève description
`Dinamis-SDK` est une bibliothèque Python qui permet de communiquer avec le prototype d'API d'infrastruture spatiale supporté par le dispositif national DINAMIS.
Le code est déposé sur le lien suivant :
https://gitlab.irstea.fr/dinamis/dinamis-sdk
%% Cell type:markdown id:f9a76706-588d-4a1c-87a5-2a236b5dfb6c tags:
### 1.2. Installation du paquet Python Dinamis-SDK
Le prototype de DINAMIS propose un catalogue d'images satellitaires dans un environnement STAC (Spatio Temporal Asset Catalogs). Nous souhaitons rechercher des images dans ce catalogue. Il faut dans un premier temps installer la bibliothèque `dinamis-sdk`.
%% Cell type:code id:dc80443b-9586-4e2b-bc33-c5f62bc72e31 tags:
``` python
#!pip install git+https://gitlab.irstea.fr/dinamis/dinamis-sdk.git
!pip install dinamis-sdk
```
%% Output
Requirement already satisfied: dinamis-sdk in /srv/conda/envs/notebook/lib/python3.10/site-packages (0.0.9)
Requirement already satisfied: requests in /srv/conda/envs/notebook/lib/python3.10/site-packages (from dinamis-sdk) (2.29.0)
Requirement already satisfied: qrcode in /srv/conda/envs/notebook/lib/python3.10/site-packages (from dinamis-sdk) (7.4.2)
Requirement already satisfied: appdirs in /srv/conda/envs/notebook/lib/python3.10/site-packages (from dinamis-sdk) (1.4.4)
Requirement already satisfied: pydantic in /srv/conda/envs/notebook/lib/python3.10/site-packages (from dinamis-sdk) (1.10.7)
Requirement already satisfied: pystac in /srv/conda/envs/notebook/lib/python3.10/site-packages (from dinamis-sdk) (1.7.3)
Requirement already satisfied: pystac-client in /srv/conda/envs/notebook/lib/python3.10/site-packages (from dinamis-sdk) (0.6.1)
Requirement already satisfied: typing-extensions>=4.2.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pydantic->dinamis-sdk) (4.5.0)
Requirement already satisfied: python-dateutil>=2.7.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pystac->dinamis-sdk) (2.8.2)
Requirement already satisfied: charset-normalizer<4,>=2 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->dinamis-sdk) (2.1.1)
Requirement already satisfied: idna<4,>=2.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->dinamis-sdk) (3.4)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->dinamis-sdk) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->dinamis-sdk) (2022.12.7)
Requirement already satisfied: pypng in /srv/conda/envs/notebook/lib/python3.10/site-packages (from qrcode->dinamis-sdk) (0.20220715.0)
Requirement already satisfied: six>=1.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from python-dateutil>=2.7.0->pystac->dinamis-sdk) (1.16.0)
%% Cell type:markdown id:88c57504-25b2-4254-a9ba-f269a3f2e4fa tags:
## 2. Découverte des images hébergées sur DINAMIS
%% Cell type:markdown id:a7204c22-be0b-4d87-aab9-2d8a16da6630 tags:
### 2.1. Accès à l'API STAC
Pour accéder au catalolgue de données, il faut importer au prélable deux bibliothèques :
- `dinamis-sdk`
- `pystac_client`
Il faut ensuite créer un objet `pystac_client.Client`.
%% Cell type:code id:803a8738-d1f2-43d3-a426-f8ef156761c6 tags:
``` python
import dinamis_sdk
import pystac_client
api = pystac_client.Client.open(
'https://stacapi-dinamis.apps.okd.crocc.meso.umontpellier.fr',
modifier=dinamis_sdk.sign_inplace,)
print(f"nom de l'api : {api.title}")
```
%% Output
nom de l'api : stac-fastapi
%% Cell type:markdown id:f54b4e1d-0976-4048-9acf-ff619eb66baa tags:
### 2.2. Recherche des collections
Dans un premier, on s'intéresse aux collections d'images disponibles dans DINAMIS. On utilise ici la méthode `get_collection()` de l'objet `api` (instance de la classe `pystac_client.client.Client`.
%% Cell type:code id:65afe82a-1601-41fe-8698-d70a898be394 tags:
``` python
collections = list(api.get_collections())
print(f"Nombre de collections: {len(collections)}")
print("Collections IDs:")
for collection in collections:
print(f"- {collection.id}")
```
%% Output
Nombre de collections: 2
Collections IDs:
- spot-6-7-drs
- super-sentinel-2-l2a
%% Cell type:markdown id:3a8cc265-79bc-412e-8a14-6e2d194cf739 tags:
### 2.3. Recherche des images
Il est désormais possible de requêter sur le catalogue afin de récupérer les items (ou images) d'intérêt. La recherche est réalisée à travers toutes les collections.
Ici, la requête repose sur deux critères :
- étendue de recherche (ou *bounding box*) et,
- une plage temporelle.
Pour chaque item, on récupère son identifiant et la collection à laquelle il appartient via les attributs `id` et `collection_id`. Ces informations permettront de retrouver aisément les images d'intérêt et de les traiter par la suite.
%% Cell type:code id:7820c73a-a0ec-4537-9893-c43060e3fed9 tags:
``` python
year = '2022'
month = '05'
bbox = [4, 42.99, 5, 44.05]
res = api.search(bbox=bbox, datetime=[f'{year}-{month}-01', f'{year}-{month}-25'])
items = res.get_all_items()
print(f"{len(items)} images répondent à la requête")
print("---")
for item in res.items():
print(f"Collection ID : {item.collection_id} ; item ID : {item.id}")
```
%% Output
12 images répondent à la requête
---
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2B_20220525-105844-883_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2B_20220522-104854-517_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2A_20220520-105851-016_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2B_20220519-103858-324_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2A_20220517-104859-725_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2A_20220514-103902-815_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2B_20220512-104852-490_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2B_20220509-103855-665_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2A_20220507-104858-836_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2A_20220504-103903-473_L2A_T31TEJ_D_15d8e3
Collection ID : super-sentinel-2-l2a ; item ID : SUPER_SENTINEL2B_20220502-104848-819_L2A_T31TEJ_D_15d8e3
Collection ID : spot-6-7-drs ; item ID : SPOT6_MS_202205011018084_SPOT6_P_202205011018084
%% Cell type:markdown id:7bbe7003-8a84-43c2-88ea-56df5a83ec28 tags:
## 3. Récupération d'une image
### 3.1. Propriétés des images
Il est possible d'accéder aux métadonnées de chaque image. Les items étant des entités geoJSON, ils peuvent être chargés sous forme de `dataframe` avec la bibliothèque `geopandas`.
%% Cell type:code id:bd88857e-fb48-49ac-a20d-aaa3ac8cd9bb tags:
``` python
import geopandas as gpd
df = gpd.GeoDataFrame.from_features(items.to_dict(), crs="epsg:4326")
df
```
%% Output
INFO:numexpr.utils:NumExpr defaulting to 8 threads.
geometry datetime
0 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-25T00:00:00Z \
1 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-22T00:00:00Z
2 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-20T00:00:00Z
3 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-19T00:00:00Z
4 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-17T00:00:00Z
5 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-14T00:00:00Z
6 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-12T00:00:00Z
7 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-09T00:00:00Z
8 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-07T00:00:00Z
9 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-04T00:00:00Z
10 POLYGON ((4.03273 43.54745, 4.03641 43.75163, ... 2022-05-02T00:00:00Z
11 POLYGON ((3.60568 44.23888, 4.33658 44.23241, ... 2022-05-01T00:00:00Z
platform instruments
0 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN] \
1 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
2 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
3 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
4 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
5 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
6 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
7 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
8 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
9 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
10 super-sentinel-2-l2a [super-sentinel-2-l2a-RGBN]
11 Spot-6 [Spot-6]
Processor_version
0 c8545bc08509415b8b505b8ae33453c889e547b8
1 c8545bc08509415b8b505b8ae33453c889e547b8
2 c8545bc08509415b8b505b8ae33453c889e547b8
3 c8545bc08509415b8b505b8ae33453c889e547b8
4 c8545bc08509415b8b505b8ae33453c889e547b8
5 c8545bc08509415b8b505b8ae33453c889e547b8
6 c8545bc08509415b8b505b8ae33453c889e547b8
7 c8545bc08509415b8b505b8ae33453c889e547b8
8 c8545bc08509415b8b505b8ae33453c889e547b8
9 c8545bc08509415b8b505b8ae33453c889e547b8
10 c8545bc08509415b8b505b8ae33453c889e547b8
11 NaN
%% Cell type:code id:e2e27816-63ec-4cd3-9a13-d17f0a6b356b tags:
``` python
#df.boundary.plot()
df.explore(column="platform",
tooltip="platform",
popup=True,
tiles="CartoDB positron",
cmap="viridis")
```
%% Output
<folium.folium.Map at 0x7f38050f8d90>
%% Cell type:markdown id:23a09599-24c9-4ce3-9ccb-42ddf10e98cc tags:
### 3.2. Sélection des items d'intérêt
On dispose d'un image Spot-6 :
- `SPOT6_MS_202205011018084_SPOT6_P_202205011018084` de la collection `spot-6-7-drs`
%% Cell type:code id:16da2d0b-3423-4366-832f-1f052b1050aa tags:
``` python
selected_spot = api.get_collection("spot-6-7-drs").get_item("SPOT6_MS_202205011018084_SPOT6_P_202205011018084")
print(selected_spot)
```
%% Output
<Item id=SPOT6_MS_202205011018084_SPOT6_P_202205011018084>
%% Cell type:markdown id:bcda050b-fd0d-4215-bcf5-d809177d9058 tags:
On cherche aussi l'image Sentinel-2 la plus proche dans le temps de l'image Spot-6.
%% Cell type:code id:f2ce29b5-232e-4301-b13f-039563e8aa13 tags:
``` python
def nearest(items, pivot):
items_s2 = [item for item in items if item.id != pivot.id]
ref = pivot.datetime
return min(items_s2, key=lambda x: abs(x.datetime - ref))
selected_s2 = nearest(items, selected_spot)
print(selected_s2)
```
%% Output
<Item id=SUPER_SENTINEL2B_20220502-104848-819_L2A_T31TEJ_D_15d8e3>
%% Cell type:markdown id:1f050d6b-8b3a-4556-a6d1-8fe8fb26f7ac tags:
### 3.3. Exploration des assets
Dans un catalogue STAC, un *item* est constitué d'*assets*. Ces derniers correspondent à un ensemble de données cohérentes entre elles. Par exemple, les *assets* d'une image SPOT-6 comprennent :
- une image multispectrale
- une image panchromatique
- des métadonnées au format Dimap
- des masques
- etc.
On récupère ici l'identifiant de l'asset et sa description avec la méthode `assets.items()` de la classe `pystac.item.Item`.
%% Cell type:code id:0f206e1e-3569-4180-91ed-00046c075b18 tags:
``` python
import rich.table
table = rich.table.Table("Asset Key", "Description")
for asset_key, asset in selected_spot.assets.items():
table.add_row(asset_key, asset.title)
table
```
%% Output
┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┓
┃ Asset Key  ┃ Description  ┃
┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━┩
│ src_xs │ │
│ src_pan │ │
│ dimap_xs │ │
│ dimap_pan │ │
│ cld_msk_vec_xs │ │
│ roi_msk_vec_xs │ │
│ cld_msk_vec_pan │ │
│ roi_msk_vec_pan │ │
│ rendered_preview │ Rendered preview │
└──────────────────┴──────────────────┘
%% Cell type:markdown id:99467500-58f6-4c60-b172-f2bbb631554b tags:
On exécute la même méthode pour l'image Sentinel-2.
%% Cell type:code id:1434ba0a-dc67-4067-ade9-ad5151684198 tags:
``` python
table = rich.table.Table("Asset Key", "Description")
for asset_key, asset in selected_s2.assets.items():
table.add_row(asset_key, asset.title)
table
```
%% Output
┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┓
┃ Asset Key  ┃ Description  ┃
┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━┩
│ clm │ │
│ img │ │
│ rendered_preview │ Rendered preview │
└──────────────────┴──────────────────┘
%% Cell type:markdown id:87bf770a-a79b-4d3f-9669-36b413e6201d tags:
### 3.4. Visualisation d'un asset
Un *asset* dispose d'un lien `href` vers la donnée. Ce `href` permet de :
- télécharger la donnée dans son ensemble,
- charger en streaming (avec le format COG) une portion de l'image
Dans un premier temps, on visualise les deux imagettes de prévisualisation. Celles-ci n'étant pas géoréférencées, il est inutile de les charger dans un client cartographique.
%% Cell type:code id:48d38a79-8a06-4e56-83fe-75fe8b0ae613 tags:
``` python
from IPython.display import Image
print("Prévisualisation RGB de l'image SPOT-6")
print("---")
Image(url=selected_spot.assets["rendered_preview"].href, width=500)
```
%% Output
Prévisualisation RGB de l'image SPOT-6
---
<IPython.core.display.Image object>
%% Cell type:code id:68559a52-11e4-4ac2-92d7-c1315b116b4d tags:
``` python
print("Prévisualisation RGB de l'image superrésolue Sentinel-2")
print("---")
Image(url=selected_s2.assets["rendered_preview"].href, width=500)
```
%% Output
Prévisualisation RGB de l'image superrésolue Sentinel-2
---
<IPython.core.display.Image object>
%% Cell type:markdown id:282bba09-acfd-48c3-8012-0c7f470c9cc6 tags:
Pour bien comprendre le contenu d'un `href`, on les affiche ici.
Les liens qui apparaissent sont non seulement "compris" dans un script Python (avec les bibliothèques qui conviennent) mais ils peuvent être également utilisés dans un logiciel SIG tel que QGis.
> Dans QGis :
> - Cliquer sur `Ajouter une couche raster`
>
> Dans la fenêtre qui apparaît :
> - Cocher le Type de source : `Protocole: HTTP(S), cloud, etc.`
> - Copier/Coller le lien href
> - Cliquer sur le bouton `Ajouter`
>
> *Dès lors la couche raster est affichée à l'écran en streaming. Elle peut être traitée avec les outils QGis (calculatrice raster, etc.).*
%% Cell type:code id:edd89178-571d-41d2-9749-b20ae8895361 tags:
``` python
print("lien href de l'image panchromatique SPOT-6")
print("---")
print(selected_spot.assets["src_pan"].href)
print("lien href de l'image Sentinel-2")
print("---")
print(selected_s2.assets["img"].href)
```
%% Output
lien href de l'image panchromatique SPOT-6
---
https://minio-api-dinamis.apps.okd.crocc.meso.umontpellier.fr/catalog/spot67/SPOT6_MS_202205011018084_SPOT6_P_202205011018084/COG_SPOT6_P_202205011018084_ORT_SPOT6_20220505_1229021hr6dt0y9l23o_1.tif?X-Amz-Security-Token=eyJhbGciOiJIUzUxMiIsInR5cCI6IkpXVCJ9.eyJhY2Nlc3NLZXkiOiIyVjlBUjhOQUlXSVRTT0lDWTJJSSIsImFjciI6IjEiLCJhbGxvd2VkLW9yaWdpbnMiOlsiKiJdLCJhdWQiOlsic2VjdXJpdHktYWRtaW4tY29uc29sZSIsImFjY291bnQiXSwiYXV0aF90aW1lIjoxNjg2NjU5ODM5LCJhenAiOiJtaW5pbyIsImVtYWlsX3ZlcmlmaWVkIjpmYWxzZSwiZXhwIjoxNjg2OTU2ODkyLCJmYW1pbHlfbmFtZSI6IktlbmppIiwiZ2l2ZW5fbmFtZSI6Ik9zZSIsImlhdCI6MTY4NjkyODA5MiwiaXNzIjoiaHR0cHM6Ly9rZXljbG9hay1kaW5hbWlzLmFwcHMub2tkLmNyb2NjLm1lc28udW1vbnRwZWxsaWVyLmZyL2F1dGgvcmVhbG1zL2RpbmFtaXMiLCJqdGkiOiJhODFlNDZiYy0xZGM4LTQ3ZGMtODNjNS00MTg5ZWI3ZGVhNGIiLCJuYW1lIjoiT3NlIEtlbmppIiwicG9saWN5IjoibXlwb2xpY3kiLCJwcmVmZXJyZWRfdXNlcm5hbWUiOiJrZW5qaS5vc2UiLCJyZWFsbV9hY2Nlc3MiOnsicm9sZXMiOlsiZGVmYXVsdC1yb2xlcy1kaW5hbWlzIiwib2ZmbGluZV9hY2Nlc3MiLCJ1bWFfYXV0aG9yaXphdGlvbiJdfSwicmVzb3VyY2VfYWNjZXNzIjp7ImFjY291bnQiOnsicm9sZXMiOlsibWFuYWdlLWFjY291bnQiLCJtYW5hZ2UtYWNjb3VudC1saW5rcyIsInZpZXctcHJvZmlsZSJdfX0sInNjb3BlIjoicHJvZmlsZSBlbWFpbCIsInNlc3Npb25fc3RhdGUiOiJkMjZjNDAzMi00M2Q3LTQyMzMtYTkwNS02ZmY3MWE2OTYxYzMiLCJzaWQiOiJkMjZjNDAzMi00M2Q3LTQyMzMtYTkwNS02ZmY3MWE2OTYxYzMiLCJzdWIiOiJhODE2ZDZhZS05NGMzLTQzMjItOTNlOS02ZGVhNzVhMDk3NjciLCJ0eXAiOiJCZWFyZXIifQ.zwEGbLIBe5rQwc1NN_KlYRLtd24NMIz-ZfMywJ4TGVkQTF2jojsgroLybBUvMMJBLorDYCOUbTTUpAtLOuFw8w&X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=2V9AR8NAIWITSOICY2II%2F20230616%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Date=20230616T150812Z&X-Amz-Expires=28800&X-Amz-SignedHeaders=host&X-Amz-Signature=a283a7fcab3ef71dd3fa3a1f39b90af363c9a012f8b81ccaea1c91ffdc8bdf0b
lien href de l'image Sentinel-2
---
https://minio-api-dinamis.apps.okd.crocc.meso.umontpellier.fr/catalog/super-sentinel-2-l2a/SUPER_SENTINEL2B_20220502-104848-819_L2A_T31TEJ_D_15d8e3/SUPER_SENTINEL2B_20220502-104848-819_L2A_T31TEJ_D_15d8e3.tif?X-Amz-Security-Token=eyJhbGciOiJIUzUxMiIsInR5cCI6IkpXVCJ9.eyJhY2Nlc3NLZXkiOiIyVjlBUjhOQUlXSVRTT0lDWTJJSSIsImFjciI6IjEiLCJhbGxvd2VkLW9yaWdpbnMiOlsiKiJdLCJhdWQiOlsic2VjdXJpdHktYWRtaW4tY29uc29sZSIsImFjY291bnQiXSwiYXV0aF90aW1lIjoxNjg2NjU5ODM5LCJhenAiOiJtaW5pbyIsImVtYWlsX3ZlcmlmaWVkIjpmYWxzZSwiZXhwIjoxNjg2OTU2ODkyLCJmYW1pbHlfbmFtZSI6IktlbmppIiwiZ2l2ZW5fbmFtZSI6Ik9zZSIsImlhdCI6MTY4NjkyODA5MiwiaXNzIjoiaHR0cHM6Ly9rZXljbG9hay1kaW5hbWlzLmFwcHMub2tkLmNyb2NjLm1lc28udW1vbnRwZWxsaWVyLmZyL2F1dGgvcmVhbG1zL2RpbmFtaXMiLCJqdGkiOiJhODFlNDZiYy0xZGM4LTQ3ZGMtODNjNS00MTg5ZWI3ZGVhNGIiLCJuYW1lIjoiT3NlIEtlbmppIiwicG9saWN5IjoibXlwb2xpY3kiLCJwcmVmZXJyZWRfdXNlcm5hbWUiOiJrZW5qaS5vc2UiLCJyZWFsbV9hY2Nlc3MiOnsicm9sZXMiOlsiZGVmYXVsdC1yb2xlcy1kaW5hbWlzIiwib2ZmbGluZV9hY2Nlc3MiLCJ1bWFfYXV0aG9yaXphdGlvbiJdfSwicmVzb3VyY2VfYWNjZXNzIjp7ImFjY291bnQiOnsicm9sZXMiOlsibWFuYWdlLWFjY291bnQiLCJtYW5hZ2UtYWNjb3VudC1saW5rcyIsInZpZXctcHJvZmlsZSJdfX0sInNjb3BlIjoicHJvZmlsZSBlbWFpbCIsInNlc3Npb25fc3RhdGUiOiJkMjZjNDAzMi00M2Q3LTQyMzMtYTkwNS02ZmY3MWE2OTYxYzMiLCJzaWQiOiJkMjZjNDAzMi00M2Q3LTQyMzMtYTkwNS02ZmY3MWE2OTYxYzMiLCJzdWIiOiJhODE2ZDZhZS05NGMzLTQzMjItOTNlOS02ZGVhNzVhMDk3NjciLCJ0eXAiOiJCZWFyZXIifQ.zwEGbLIBe5rQwc1NN_KlYRLtd24NMIz-ZfMywJ4TGVkQTF2jojsgroLybBUvMMJBLorDYCOUbTTUpAtLOuFw8w&X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=2V9AR8NAIWITSOICY2II%2F20230616%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Date=20230616T150812Z&X-Amz-Expires=28800&X-Amz-SignedHeaders=host&X-Amz-Signature=98316a8f0267eda512ad72c34d0c6db518ff5b8eed951f62528d6a678b9da042
%% Cell type:markdown id:b39e1bd8-184f-4bab-8093-57fce558c33e tags:
## 4. Calcul du NDVI
Pour finir, on souhaite calculer l'indice NDVI sur une portion de l'image SPOT-6. Il faut donc récupérer une portion (définie par une emprise, ou *bounding box*) de l'image multipectrale, et réaliser un calcul entre bandes.
%% Cell type:markdown id:dbbac706-b83e-4cac-b4fe-164fe42ee360 tags:
### 4.1 Préparation de l'emprise de découpage
L'emprise de découpage est un rectangle définie par les coordonnées extrémales suivantes :
- X min : 775000,
- Y min : 6328000,
- X max : 780000,
- Y max : 6333000
Celles-ci sont exprimées en mètres dans le système de coordonnées RGF-93 / Lambert-93 (EPSG:2154).
On transforme tout d'abord ces coordonnées en objet geoJSON. Ce type est attendu comme argument lors du découpage.
%% Cell type:code id:9c908eac-6452-4630-8293-0849b76d0f15 tags:
``` python
from shapely.geometry import box
from fiona.crs import from_epsg
bbox = box(775000, 6328000, 780000, 6333000)
geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=from_epsg(2154))
def getFeatures(gdf):
"""Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""
import json
return [json.loads(gdf.to_json())['features'][0]['geometry']]
coords = getFeatures(geo)
coords
```
%% Output
[{'type': 'Polygon',
'coordinates': [[[780000.0, 6328000.0],
[780000.0, 6333000.0],
[775000.0, 6333000.0],
[775000.0, 6328000.0],
[780000.0, 6328000.0]]]}]
%% Cell type:markdown id:9dddf806-129b-40c5-96ed-3e0bc881540d tags:
### 4.2. Chargement d'une portion d'image et calcul du NDVI
On utilise la bibliothèque `rasterio` pour lire l'image multispectrale (référencée par le `href`). La méthode `mask` associée permet de requêter uniquement sur la portion d'image définie par la bounding box (variable `coords`).
On procède ensuite au calcul du NDVI (*Normalized Difference Vegetation Index*), en suivant l'expression suivante :
$\dfrac{(pIR - rouge)}{(pIR + rouge)}$
%% Cell type:code id:d52ec352-388c-457f-b5fc-e7b50fa2c9e0 tags:
``` python
import rasterio as rio
from rasterio.mask import mask
import numpy as np
with rio.open(selected_spot.assets["src_xs"].href) as xs:
#print(xs.crs)
clipXs, out_transform = mask(dataset=xs, shapes=coords, crop=True)
red = clipXs[0]
nir = clipXs[3]
ndvi = np.zeros(red.shape, dtype=rio.float32)
ndvi = (nir-red)/(nir+red+0.0001)
```
%% Cell type:markdown id:e3783c3b-1331-4140-ab9f-59a38bec8f11 tags:
### 4.3. Affichage de l'image NDVI
On propose ici une autre manière de visualiser un raster, en passant par le module `pyplot` de `matplotlib`.
%% Cell type:code id:c4275daa-6a48-40e4-8d1c-f1ee9285ede2 tags:
``` python
import matplotlib.pyplot as plt
plt.imshow(ndvi, cmap="RdYlGn")
plt.colorbar()
plt.show()
```
%% Cell type:markdown id:3ae302df-df13-47be-b26f-223c2938ab3a tags:
### 4.4. Utilisation des widgets Jupyter
Il est possible d'ajouter des éléments (*slider*, *progress bar*, *checkbox*, *radio buttons*, etc.) pour interagir avec la visualisation des données. Pour ce faire, charger la librairie `ipywidgets`.
%% Cell type:code id:1bca7f81-8cb7-4f00-9989-824ede9b2620 tags:
``` python
import ipywidgets as widgets
from ipywidgets import interact, interact_manual
def threshold(seuil):
seuil = np.where(ndvi>seuil, 1, 0)
plt.imshow(seuil, cmap="viridis", interpolation='nearest')#, cmap=plt.cm.gray)
plt.colorbar()
plt.show()
interact(threshold, seuil = widgets.FloatSlider(min=-1, max=1, step=0.001))
```
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment