Extraction de valeurs à une latitude et une longitude spécifiques à partir des données d'andain MODIS

J'essaie de déterminer la quantité de vapeur d'eau précipitable (PWV), d'ozone et d'aérosols en fonction du temps sur des points spécifiques de la Terre, à savoir nos observatoires astronomiques. Pour ce faire, j'ai déjà du code Python modapsclientqui télécharge les produits MODIS Aqua et Terra MYDATML2 et MODATML2 deux fois par jour qui couvrent la latitude et la longitude spécifiques qui m'intéressent.

Ce que je ne sais pas, c'est comment extraire les quantités spécifiques que je veux, telles que l'heure à laquelle les données MODIS ont été prises et PWV pour la position particulière de latitude et de longitude de mon observatoire pour en faire une série temporelle de valeurs. Les produits MYDATML2 semblent contenir des grilles de latitude et de longitude 2D Cell_Along_Swath_5kmet Cell_Across_Swath_5kmdonc je suppose que cela en fait des données en andains plutôt que des données de mosaïque ou de grille? Les quantités que je veux telles que Precipitable_Water_Infrared_ClearSkysemblent également être contre Cell_Along_Swath_5kmet mais Cell_Across_Swath_5kmje ne sais pas comment obtenir cette valeur PWV au lat spécifique, longtemps je suis intéressé. Aide s'il vous plaît?

[EDIT 1 - J'ai changé la recherche de coordonnées de pixels]

En utilisant cet exemple de MODATML que vous avez fourni et en utilisant la bibliothèque gdal. Ouvrons le hdf avec gdal:

import gdal
dataset = gdal.Open(r"E:\modis\MODATML2.A2018182.0800.061.2018182195418.hdf")

Ensuite, nous voulons voir comment les sous-jeux de données sont nommés afin d'importer correctement ceux dont nous avons besoin:

datasets_meta = dataset.GetMetadata("SUBDATASETS")

Cela renvoie un dictionnaire:

datasets_meta
>>>{'SUBDATASET_1_NAME': 'HDF4_EOS:EOS_SWATH:"E:\\modis\\MODATML2.A2018182.0800.061.2018182195418.hdf":atml2:Cloud_Optical_Thickness', 
'SUBDATASET_1_DESC': '[406x271] Cloud_Optical_Thickness atml2 (16-bit integer)',
'SUBDATASET_2_NAME':'HDF4_EOS:EOS_SWATH:"E:\\modis\\MODATML2.A2018182.0800.061.2018182195418.hdf":atml2:Cloud_Effective_Radius',
'SUBDATASET_2_DESC': '[406x271] Cloud_Effective_Radius atml2 (16-bit integer)',
[....]}

Disons que nous voulons obtenir la première variable, l'épaisseur optique du nuage, nous pouvons accéder à son nom par:

datasets_meta['SUBDATASET_1_NAME']
>>>'HDF4_EOS:EOS_SWATH:"E:\\modis\\MODATML2.A2018182.0800.061.2018182195418.hdf":atml2:Cloud_Optical_Thickness'

Maintenant, nous pouvons charger la variable en mémoire en appelant à nouveau la méthode .Open ():

Cloud_opt_th = gdal.Open(datasets_meta['SUBDATASET_1_NAME'])

Par exemple, vous pouvez accéder à Precipitable_Water_Infrared_ClearSky qui vous intéresse en fournissant «SUBDATASET_20_NAME». Jetez un œil au dictionnaire datasets_meta.

Cependant, la variable extraite n'a pas de géoprojection (var.GetGeoprojection ()) comme vous vous en doutez pour d'autres types de fichiers tels que GeoTiff. Vous pouvez charger la variable sous forme de tableau numpy et tracer la variable 2d sans projection:

Cloud_opt_th_array = Cloud_opt_th.ReadAsArray()
import matplotlib.pyplot as plt
plt.imshow(Cloud_opt_th_array)

Maintenant, puisqu'il n'y a pas de géoprojection, nous allons examiner les métadonnées de la variable:

Cloud_opt_th_meta = Cloud_opt_th.GetMetadata()

Ceci est un autre dictionnaire qui comprend toutes les informations dont vous avez besoin, y compris une longue description du sous-échantillonnage (j'ai remarqué que cela n'est fourni qu'avec le premier sous-ensemble de données), qui comprend l'explication de ces Cell_Along_Swath:

Cloud_opt_th_meta['1_km_to_5_km_subsampling_description']
>>>'Each value in this dataset does not represent an average of properties over a 5 x 5 km grid box, but rather a single sample from within each 5 km box. Normally, pixels in across-track rows 4 and 9 (counting in the direction of increasing scan number) out of every set of 10 rows are used for subsampling the 1 km retrievals to a 5 km resolution. If the array contents are determined to be all fill values after selecting the default pixel subset (e.g., from failed detectors), a different pair of pixel rows is used to perform the subsampling. Note that 5 km data sets are centered on rows 3 and 8; the default sampling choice of 4 and 9 is for better data quality and avoidance of dead detectors on Aqua. The row pair used for the 1 km sample is always given by the first number and last digit of the second number of the attribute Cell_Along_Swath_Sampling. The attribute Cell_Across_Swath_Sampling indicates that columns 3 and 8 are used, as they always are, for across-track sampling. Again these values are to be interpreted counting in the direction of the scan, from 1 through 10 inclusively. For example, if the value of attribute Cell_Along_Swath_Sampling is 3, 2028, 5, then the third and eighth pixel rows were used for subsampling. A value of 4, 2029, 5 indicates that the default fourth and ninth rows pair was used.'

Je pense que cela signifie que sur la base de ces pixels de 1 km, le 5 km a été construit en prenant exactement les valeurs des pixels à une certaine position dans le tableau de détection 5x5 (la position est indiquée dans les métadonnées, je pense que c'est une chose instrumentale pour réduire les défauts).

Quoi qu'il en soit, à ce stade, nous avons un tableau de cellules de 1x1 km (voir la description du sous-échantillonnage ci-dessus, pas sûr de la science derrière). Pour obtenir les coordonnées de chaque pixel centroïde, nous devons charger les sous-jeux de latitude et de longitude.

Latitude = gdal.Open(datasets_meta['SUBDATASET_66_NAME']).ReadAsArray()
Longitude = gdal.Open(datasets_meta['SUBDATASET_67_NAME']).ReadAsArray()

Par exemple,

Longitude
>>> array([[-133.92064, -134.1386 , -134.3485 , ..., -154.79303, -154.9963 ,
    -155.20723],
   [-133.9295 , -134.14743, -134.3573 , ..., -154.8107 , -155.01431,
    -155.2256 ],
   [-133.93665, -134.1547 , -134.36465, ..., -154.81773, -155.02109,
    -155.23212],
   ...,
   [-136.54477, -136.80055, -137.04684, ..., -160.59378, -160.82101,
    -161.05663],
   [-136.54944, -136.80536, -137.05179, ..., -160.59897, -160.8257 ,
    -161.06076],
   [-136.55438, -136.81052, -137.05714, ..., -160.6279 , -160.85527,
    -161.09099]], dtype=float32)        

Vous pouvez remarquer que les coordonnées Latitude et Longitude sont différentes pour chaque pixel.

Supposons que votre observatoire se trouve aux coordonnées lat_obs, long_obs, puis minimisez la différence de coordonnées x, y:

coordinates = np.unravel_index((np.abs(Latitude - lat_obs) + np.abs(Longitude - long_obs)).argmin(), Latitude.shape)

et extrayez votre valeur

Cloud_opt_th_array[coordinates]