Comment lire un fichier MODIS HDF en utilisant Python et pyhdf ?

Published: 15 octobre 2024

Tags: Python; MODIS; HDF;

DMCA.com Protection Status

Introduction

MODIS est un instrument clé des satellites Terra et Aqua de la NASA qui mesure les dynamiques globales à grande échelle, y compris la couverture nuageuse, le rayonnement et d'autres données atmosphériques. Ces données sont stockées au format HDF4, qui est légèrement différent du format HDF5 plus courant.

Ce tutoriel vous guidera à travers le processus d'utilisation de la bibliothèque pyhdf en Python pour lire un fichier HDF (Hierarchical Data Format) MODIS (Moderate Resolution Imaging Spectroradiometer). Plus précisément, nous travaillerons avec un fichier MODIS de niveau 2, qui contient des ensembles de données scientifiques (SDS) liés aux propriétés des nuages.

Exigences

Assurez-vous d'avoir les packages suivants installés dans votre environnement Python :

  • pyhdf : pour lire les fichiers HDF4 (Installez via pip install pyhdf ou conda install -c anaconda pyhdf)
  • numpy : pour la manipulation de tableaux (Installez via pip install numpy)

Téléchargez le fichier HDF MODIS avec lequel nous allons travailler à partir de ce lien.

De plus, vous pouvez télécharger une version Jupyter notebook de ce guide ici.

Lecture d'un fichier HDF4 MODIS

Une fois le fichier téléchargé, commencez par importer les modules nécessaires et ouvrir le fichier HDF :

1
2
3
4
5
from pyhdf.SD import SD, SDC
import numpy as np

file_name = './inputs/MYD06_L2.A2007219.2010.006.2014053202546.hdf'
file = SD(file_name, SDC.READ)

La classe SD est utilisée pour ouvrir le fichier en mode lecture. Le drapeau SDC.READ l'ouvre en mode lecture seule.

Pour vérifier les informations générales sur le fichier, telles que le nombre d'ensembles de données scientifiques (SDS) qu'il contient, utilisez :

1
file.info()

Vous devriez voir que le fichier contient 127 SDS. Chaque SDS représente une mesure scientifique différente, telle que la température du sommet des nuages ou le masque nuageux.

Impression des ensembles de données scientifiques (SDS)

Pour lister tous les noms des SDS, nous pouvons parcourir la méthode datasets() du fichier :

1
2
3
4
datasets_dic = file.datasets()

for idx, sds in enumerate(datasets_dic.keys()):
    print(idx, sds)

La sortie affichera tous les SDS dans le fichier, tels que :

1
2
3
4
5
6
0 Cloud_Top_Height_Nadir_Night
1 Single_Scatter_Albedo_Ice
2 Radiance_Variance
3 Brightness_Temperature
...
126 Cloud_Effective_Radius_Uncertainty_1621

Accès aux données SDS spécifiques

Une fois que nous avons identifié le SDS avec lequel nous voulons travailler, nous pouvons extraire ses données. Par exemple, pour récupérer les données du SDS Cloud Top Temperature at 1km, utilisez :

1
2
3
sds_obj = file.select('cloud_top_temperature_1km')  # Sélectionner le SDS
data = sds_obj.get()  # Obtenir les données du SDS
print(data)

Cela affichera les données brutes sous forme de tableau :

1
2
3
4
[[11758 11565 10858 ... 12100 12299 11902]
 [11758 11565 10667 ... 12100 12299 12100]
 [11758 11565 10858 ... 11902 12299 12100]
 ...]

Interprétation des attributs SDS

Chaque SDS contient des métadonnées, telles que des unités et des facteurs d'échelle, qui aident à interpréter les données brutes. Pour accéder aux attributs du SDS cloud_top_temperature_1km, utilisez :

1
2
import pprint
pprint.pprint(sds_obj.attributes())

La sortie ressemblera à ceci :

1
2
3
4
5
6
{'_FillValue': -999,
 'add_offset': -15000.0,
 'long_name': 'Cloud Top Temperature at 1-km resolution from LEOCAT, '
 'scale_factor': 0.009999999776482582,
 'units': 'K',
 'valid_range': [0, 20000]}

Les attributs add_offset et scale_factor sont essentiels pour convertir les données brutes en valeurs significatives. Pour les appliquer aux données :

1
2
3
4
5
6
add_offset = sds_obj.attributes()['add_offset']
scale_factor = sds_obj.attributes()['scale_factor']

# Reformater les données
data = (data - add_offset) * scale_factor
print(data)

Cela convertira les valeurs entières brutes en températures réelles en Kelvin (K).

Extraction de données complexes (ex. : masque nuageux)

Certaines SDS stockent plusieurs informations dans un seul tableau de données. Par exemple, le SDS Cloud_Mask_1km stocke plusieurs indicateurs dans son format compact par bits.

1
2
3
sds_obj = file.select('Cloud_Mask_1km')
data = sds_obj.get()
print(data.shape)  # (2030, 1354, 2)

Pour extraire des indicateurs spécifiques, tels que l'indicateur d'état du masque nuageux, l'indicateur de nébulosité du masque nuageux et l'indicateur jour/nuit, utilisez une fonction de stripping des bits :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def bits_stripping(bit_start, bit_count, value):
    bitmask = pow(2, bit_start + bit_count) - 1
    return np.right_shift(np.bitwise_and(value, bitmask), bit_start)

i, j = 576, 236  # Pixel aléatoire

status_flag = bits_stripping(0, 1, data[i, j, 0]) 
day_flag = bits_stripping(3, 1, data[i, j, 0]) 
cloud_mask_flag = bits_stripping(1, 2, data[i, j, 0])

print('Masque nuageux déterminé:', status_flag)
print('Indicateur jour/nuit:', day_flag)
print('Masque nuageux:', cloud_mask_flag)

Cela décodera les informations compactées dans l'octet.

Code Python complet (Python 3)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#!/usr/bin/env python

from pyhdf.SD import SD, SDC
import numpy as np
import pprint

file_name = 'MYD06_L2.A2007219.2010.006.2014053202546.hdf'
file = SD(file_name, SDC.READ)

# Imprimer des informations de base
print(file.info())

# Lister les noms des SDS
datasets_dic = file.datasets()
for idx, sds in enumerate(datasets_dic.keys()):
    print(idx, sds)

# Accéder et traiter le SDS de la température du sommet des nuages
sds_obj = file.select('cloud_top_temperature_1km')
data = sds_obj.get()

# Extraire et appliquer le scale_factor et l'add_offset
attrs = sds_obj.attributes()
add_offset = attrs['add_offset']
scale_factor = attrs['scale_factor']
data = (data - add_offset) * scale_factor
print(data)

# Exemple de masque nuageux
sds_obj = file.select('Cloud_Mask_1km')
data = sds_obj.get()
print(data.shape)

# Extraire les informations compactées par bits
def bits_stripping(bit_start, bit_count, value):
    bitmask = pow(2, bit_start + bit_count) - 1
    return np.right_shift(np.bitwise_and(value, bitmask), bit_start)

i, j = 576, 236
status_flag = bits_stripping(0, 1, data[i, j, 0])
day_flag = bits_stripping(3, 1, data[i, j, 0])
cloud_mask_flag = bits_stripping(1, 2, data[i, j, 0])

print('Masque nuageux déterminé:', status_flag)
print('Indicateur jour/nuit:', day_flag)
print('Masque nuageux:', cloud_mask_flag)

Compatibilité Python 2.7

Si vous utilisez Python 2.7, le code est presque identique. Cependant, vous devrez ajuster les instructions d'impression et l'itération sur les dictionnaires :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#!/usr/bin/env python

# Importation des modules nécessaires de pyhdf pour lire les fichiers HDF4
from pyhdf.SD import SD, SDC

# Spécifiez le nom du fichier HDF à lire (assurez-vous que le chemin du fichier est correct)
file_name = 'MYD06_L2.A2007219.2010.006.2014053202546.hdf'

# Ouvrir le fichier HDF en lecture (SDC.READ ouvre le fichier en mode lecture seule)
file = SD(file_name, SDC.READ)

# Imprimer les informations générales du fichier, y compris le nombre de jeux de données, les dimensions, etc.
print file.info()

# --------------------------------------------------------------------------------

--
# Lister les jeux de données scientifiques (SDS) présents dans le fichier HDF
datasets_dic = file.datasets()
for idx, sds in datasets_dic.items():
    print idx, sds

# Sélectionner et traiter un SDS spécifique : température du sommet des nuages
sds_obj = file.select('cloud_top_temperature_1km')  # Sélectionner le SDS
data = sds_obj.get()  # Obtenir les données du SDS

# Accéder aux attributs du SDS
attrs = sds_obj.attributes()
add_offset = attrs['add_offset']
scale_factor = attrs['scale_factor']

# Reformater les données
data = (data - add_offset) * scale_factor
print data

# Exemple d'accès au masque nuageux
sds_obj = file.select('Cloud_Mask_1km')
data = sds_obj.get()
print data.shape

# Fonction pour extraire les indicateurs à partir d'un tableau de données compact
def bits_stripping(bit_start, bit_count, value):
    bitmask = pow(2, bit_start + bit_count) - 1
    return np.right_shift(np.bitwise_and(value, bitmask), bit_start)

# Exemple d'indice pour accéder à un pixel spécifique
i, j = 576, 236
status_flag = bits_stripping(0, 1, data[i, j, 0])
day_flag = bits_stripping(3, 1, data[i, j, 0])
cloud_mask_flag = bits_stripping(1, 2, data[i, j, 0])

# Imprimer les indicateurs extraits
print 'Masque nuageux déterminé:', status_flag
print 'Indicateur jour/nuit:', day_flag
print 'Masque nuageux:', cloud_mask_flag

Références

Liens Source
MODIS C6 cloud optical products user guide NASA
pyhdf sourceforge