Opening EOS netCDF4/HDF5 file with correct format using xarray?

Question

I am using xarray (version 0.12.2) to work with the MODIS/Terra+Aqua MAIAC Land Aerosol Optical Depth dataset (MCD19A2.006). xarray successfully opens that dataset with the correct variables and dimensions but the dataset has projection metadata that is currently not being read. Has anyone else encountered this problem or have solutions to parsing the metadata to construct the correct coordinates?

Here is an example:

import xarray as xr

ds = xr.open_dataset('MCD19A2.A2000057.h09v07.006.2018013034454.hdf')

print(ds)

Which outputs the following dataset information:

<xarray.Dataset>
Dimensions:            (Orbits:grid1km: 1, Orbits:grid5km: 1, XDim:grid1km: 1200, XDim:grid5km: 240, YDim:grid1km: 1200, YDim:grid5km: 240)
Dimensions without coordinates: Orbits:grid1km, Orbits:grid5km, XDim:grid1km, XDim:grid5km, YDim:grid1km, YDim:grid5km
Data variables:
    Optical_Depth_047  (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Optical_Depth_055  (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_Uncertainty    (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    FineModeFraction   (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Column_WV          (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_QA             (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_MODEL          (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Injection_Height   (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    cosSZA             (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    cosVZA             (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    RelAZ              (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    Scattering_Angle   (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    Glint_Angle        (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
Attributes:
    HDFEOSVersion:                     HDFEOS_V2.19
    StructMetadata.0:                  GROUP=SwathStructure\nEND_GROUP=SwathS...
    Orbit_amount:                      1
    Orbit_time_stamp:                  20000571645T  
    CoreMetadata.0:                    \nGROUP                  = INVENTORYME...
    ArchiveMetadata.0:                 \nGROUP                  = ARCHIVEDMET...
    identifier_product_doi:            10.5067/MODIS/MCD19A2.006
    identifier_product_doi_authority:  http://dx.doi.org

You can see that there are no coordinates associated with dataset but in reality it does have geographic coordinates, it is just buried deep in metadata and needs to be extracted. My next step is to read the projection information from the StructMetadata.0 attribute and construct the coordinates from there. But I wonder if I am missing something with xarray or if someone else has code to where it would extract the geographic information automatically.

Answer 1

I ended up writing a function to parse the structured metadata attribute and add coordinates to a dataset. Here is my code to read in the dataset, create the coordinates, and apply it to the dataset:

import xarray as xr
import pyproj
import numpy as np
from collections import OrderedDict

def parse_hdfeos_metadata(string):
  out = OrderedDict()
  lines = [i.replace('\t','') for i in string.split('\n')]
  i = -1
  while i<(len(lines))-1:
      i+=1
      line = lines[i]
      if "=" in line:
          key,value = line.split('=')
          if key in ['GROUP','OBJECT']:
              endIdx = lines.index('END_{}={}'.format(key,value))
              out[value] = parse_hdfeos_metadata("\n".join(lines[i+1:endIdx]))
              i = endIdx
          else:
              if ('END_GROUP' not in key) and ('END_OBJECT' not in key):
                   try:
                       out[key] = eval(value)
                   except NameError:
                       out[key] = str(value)
  return out

def construct_coords(ds,grid='GRID_1'):
    metadata = parse_hdfeos_metadata(ds.attrs['StructMetadata.0'])

    gridInfo = metadata['GridStructure'][grid]

    gridName = gridInfo['GridName']

    x1,y1 = gridInfo['UpperLeftPointMtrs']
    x2,y2 = gridInfo['LowerRightMtrs']
    yRes = (y1-y2)/gridInfo['YDim']
    xRes = (x1-x2)/gridInfo['XDim']

    x = np.arange(x2,x1,xRes)
    y = np.arange(y2,y1,yRes)[::-1]

    xx,yy = np.meshgrid(x,y)

    # hard coded projection name information...
    if 'soid' in gridInfo['Projection'].lower():
        pp = 'sinu'
    else:
        pp = gridInfo['Projection'].lower()

    projStr = "+proj={} +lon_0=0 +x_0=0 +y_0=0 +a={} +units=m +no_defs".format(
      pp,gridInfo['ProjParams'][0])

    proj = pyproj.Proj(projStr)
    gcs = proj.to_latlong()

    lon,lat = pyproj.transform(proj,gcs,xx,yy)

    ydim,xdim = 'YDim:{}'.format(gridName),'XDim:{}'.format(gridName)
    yCoordName = 'Latitude:{}'.format(gridName)
    xCoordName = 'Longitude:{}'.format(gridName)

    out = ds.copy()

    out[yCoordName], out[xCoordName] = ((ydim,xdim),lat),((ydim,xdim),lon)

    return out.set_coords([yCoordName,xCoordName])


ds = xr.open_dataset('MCD19A2.A2000057.h09v07.006.2018013034454.hdf')
geo_ds = construct_coords(ds)

print(geo_ds)

With the following output that has coordinates:

<xarray.Dataset>
Dimensions:            (Orbits:grid1km: 4, Orbits:grid5km: 4, XDim:grid1km: 1200, XDim:grid5km: 240, YDim:grid1km: 1200, YDim:grid5km: 240)
Coordinates:
    Latitude:grid1km   (YDim:grid1km, XDim:grid1km) float64 19.99 19.99 ... 10.0
    Longitude:grid1km  (YDim:grid1km, XDim:grid1km) float64 85.13 ... 71.09
Dimensions without coordinates: Orbits:grid1km, Orbits:grid5km, XDim:grid1km, XDim:grid5km, YDim:grid1km, YDim:grid5km
Data variables:
    Optical_Depth_047  (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Optical_Depth_055  (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_Uncertainty    (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    FineModeFraction   (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Column_WV          (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_QA             (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    AOD_MODEL          (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    Injection_Height   (Orbits:grid1km, YDim:grid1km, XDim:grid1km) float32 ...
    cosSZA             (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    cosVZA             (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    RelAZ              (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    Scattering_Angle   (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
    Glint_Angle        (Orbits:grid5km, YDim:grid5km, XDim:grid5km) float32 ...
Attributes:
    HDFEOSVersion:                     HDFEOS_V2.19
    StructMetadata.0:                  GROUP=SwathStructure\nEND_GROUP=SwathS...
    Orbit_amount:                      4
    Orbit_time_stamp:                  20190220450T  20190220630T  2019022075...
    CoreMetadata.0:                    \nGROUP                  = INVENTORYME...
    ArchiveMetadata.0:                 \nGROUP                  = ARCHIVEDMET...
    identifier_product_doi:            10.5067/MODIS/MCD19A2.006
    identifier_product_doi_authority:  http://dx.doi.org