Extracting values at specific latitude,longitude from MODIS swath data

Question

I am trying to determine the amount of precipitable water vapor (PWV), ozone and aerosols as a function of time over specific spots on the Earth, namely our astronomical observatories. To do this, I've already got some Python code using modapsclient which will download the twice daily MODIS Aqua and Terra MYDATML2 and MODATML2 products that cover the specific latitude and longitude I am interested in.

What I'm not sure about is how to extract the specific quantities I want such as the time the MODIS data were taken and PWV for the particular latitude and longitude position of my observatory to make them into a time series of values. The MYDATML2 products seem to contain 2D latitude and longitude grids of Cell_Along_Swath_5km and Cell_Across_Swath_5km so I guess this makes it swath data as opposed to tile or grid data ? The quantities I want such as Precipitable_Water_Infrared_ClearSky also seem to be against the Cell_Along_Swath_5km and Cell_Across_Swath_5km butI'm not sure how to get that PWV value at the specific lat,long I am interested in. Help please ?

Answer 1

[EDIT 1 - I changed pixel coordinate search]

Using this sample of MODATML that you provided and using gdal library. Let's open the hdf with gdal:

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

Then we want to view how the subdatasets are named in order to correctly import the ones we need:

datasets_meta = dataset.GetMetadata("SUBDATASETS")

This returns a dictionary:

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)',
[....]}

Let's say we want to get the first variable, the cloud optical thickness, we can access its name by:

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

Now we can load the variable in memory calling again .Open() method:

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

For example, you can access Precipitable_Water_Infrared_ClearSky you are interested in by providing 'SUBDATASET_20_NAME'. Just have a look at datasets_meta dictionary.

However, the variable extracted does not have a geoprojection (var.GetGeoprojection()) as you would expect from other file types such as GeoTiff. You can load the variable as a numpy array and plot the 2d variable without projection:

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

Now, since there is no geoprojection, we will look into metadata of the variable:

Cloud_opt_th_meta = Cloud_opt_th.GetMetadata()

This is another dictionary That includes all the information you need, including a long description of the subsampling(I noticed this is provided only with the first subdataset), that includes the explanation of these 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.'

I think this means that based on these 1km pixels the 5km was built taking exactly the pixels values at a certain position in the 5x5 sensing array (the position is indicated in the metadata, I think this is an instrument thing to reduce faults).

Anyhow, at this point we have an array of cells 1x1 km (see description of subsampling above, not sure about the science behind it). To get the coordinates of each pixel centroid, we need to load the latitude and longitude subdatasets.

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

For example,

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)        

You may notice that Latitude and Longitude coordinates are different for each pixel.

Say your observatory is located at coordinates lat_obs, long_obs, than you minimise the x,y coordinate difference:

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

and extract your value

Cloud_opt_th_array[coordinates]