[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:
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:
'SUBDATASET_1_DESC': '[406x271] Cloud_Optical_Thickness atml2 (16-bit integer)',
'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:
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
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:
>>>'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()
>>> array([[-133.92064, -134.1386 , -134.3485 , ..., -154.79303, -154.9963 ,
[-133.9295 , -134.14743, -134.3573 , ..., -154.8107 , -155.01431,
[-133.93665, -134.1547 , -134.36465, ..., -154.81773, -155.02109,
[-136.54477, -136.80055, -137.04684, ..., -160.59378, -160.82101,
[-136.54944, -136.80536, -137.05179, ..., -160.59897, -160.8257 ,
[-136.55438, -136.81052, -137.05714, ..., -160.6279 , -160.85527,
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