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.

Colabean · Accepted Answer · 2020-03-09 15:31:09Z

To open the data with the projection information you need to open the sub-datasets individually.

I will use a MODIS dataset I have to hand as an example, MOD11A1, but it will be the same for yours. You can get the filename of the subdatasets using rasterio for example:

import rasterio
filename = '/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf'
with rasterio.open(filename) as src:
    subdatasets = src.subdatasets

You could use gdal rather than rasterio:

import gdal
g = gdal.Open(filename)
subdatasets = g.GetSubDatasets()

In this example, subdatasets looks like:

print(subdatasets)
['HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:Emis_32', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:Clear_day_cov', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:Clear_night_cov', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:QC_Day', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:Day_view_time', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:Day_view_angl', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:LST_Night_1km', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:QC_Night', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:Night_view_time', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:Night_view_angl', 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:Emis_31']

Opening one of these subdatasets as an xarray will preserve the projection information:

import xarray as xr
fname = 'HDF4_EOS:EOS_GRID:/data/MOD11A1.A2019225.h17v03.006.2019226085002.hdf:MODIS_Grid_Daily_1km_LST:LST_Day_1km'
myDataset = xr.open_rasterio(fname)

And I have an xarray with projection information:

print(myDataset)
<xarray.DataArray (band: 1, y: 1200, x: 1200)>
[1440000 values with dtype=uint16]
Coordinates:
  * band     (band) int64 1
  * y        (y) float64 6.671e+06 6.67e+06 6.669e+06 ... 5.561e+06 5.56e+06
  * x        (x) float64 -1.111e+06 -1.111e+06 -1.11e+06 ... -1.39e+03 -463.3
Attributes:
    transform:     (926.6254331391667, 0.0, -1111950.519767, 0.0, -926.625433...
    crs:           +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=637100...
    res:           (926.6254331391667, 926.6254331383334)
    is_tiled:      0
    nodatavals:    (0.0,)
    scales:        (0.02,)
    offsets:       (0.0,)
    descriptions:  ('Daily daytime 1km grid Land-surface Temperature',)
    units:         ('K',)

If you need all of the sub-datasets it is necessary to loop though each sub-product and then add them to an xarray dataset.

Kel Markert · Accepted Answer · 2019-07-19 18:25:32Z

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

I have been trying with some success to adjust this code to read and display MODIS data as well. Just wondering if you ended up with a better solution. — mmann1123, Jan 29 at 20:18
I have not come up with another solution, this has been sufficient for my purposes. What exactly are you trying to do with the MODIS data (and which dataset)? Happy to help adjust the code if needed. — Kel Markert, Jan 30 at 15:15

Stack Exchange Network

current community

your communities

more stack exchange communities

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

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged python hdf5 xarray or ask your own question.

Hot Network Questions

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

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged python hdf5 xarray or ask your own question.

Related

Hot Network Questions