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.
Table des matières
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
ouconda 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 |