xarray wrapper for MODIS

xr_modis.py

__all__ = ['open_modis']
__version__ = '0.1.0'


def parse_hdfeos_metadata(string):
    """
    Borrowed from ARSET training materials
    """
    out = dict()
    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 open_modis(path, keys=None, verbose=0):
    from pyhdf import SD
    import xarray as xr
    import pandas as pd
    import numpy as np

    sds = SD.SD(path)
    attrs = sds.attributes()
    dadefs = sds.datasets()
    meta = parse_hdfeos_metadata(attrs['StructMetadata.0'])
    earth_radius = meta['GridStructure']['GRID_1']['ProjParams'][0]
    attrs['crs_proj4'] = f"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a={earth_radius} +units=m +no_defs"
    attrs.pop('StructMetadata.0', '')
    attrs.pop('CoreMetadata.0', '')
    attrs.pop('ArchiveMetadata.0', '')
    ds = xr.Dataset(attrs=attrs)
    if keys is None:
        keys = list(dadefs)
    for key, (dims, shape, dtype, index) in dadefs.items():
        if key not in keys:
            continue
        sdv = sds.select(key)
        sdvattrs = sdv.attributes()
        sdval = sdv.get() * sdvattrs.get('scale_factor', 1) + sdvattrs.get('add_offset', 0)
        if verbose > 0:
            print(key)
        ds[key] = xr.DataArray(sdval, dims=dims, attrs=sdvattrs)

    for gridkey, griddef in meta['GridStructure'].items():
        gdnam = griddef['GridName']
        xkey = f'XDim:{gdnam}'
        ykey = f'YDim:{gdnam}'
        x = np.linspace(griddef['UpperLeftPointMtrs'][0], griddef['LowerRightMtrs'][0], griddef['XDim'])
        y = np.linspace(griddef['LowerRightMtrs'][1], griddef['UpperLeftPointMtrs'][1], griddef['YDim'])
        ds.coords[xkey] = x
        ds.coords[ykey] = y
        ds.coords[f'Orbits:{gdnam}'] = pd.to_datetime(
            [t[:-1] for t in attrs['Orbit_time_stamp'].split()], format='%Y%j%H%M'
        )
    return ds