Module pyl4c.lib.modis

MODIS sinusoidal projection forward and backward coordinate transformations, courtesy of Giglio et al. (2018), Collection 6 MODIS Burned Area Product User's Guide, Version 1, Appendix B:

<https://modis-fire.umd.edu/files/MODIS_C6_BA_User_Guide_1.2.pdf>

Functions

def dec2bin_unpack(x)
Expand source code
def dec2bin_unpack(x):
    '''
    Unpacks an arbitrary decimal NumPy array into a binary representation
    along a new axis. Assumes decimal digits are on the interval [0, 255],
    i.e., that only 8-bit representations are needed.

    Parameters
    ----------
    x : numpy.ndarray

    Returns
    -------
    numpy.ndarray
    '''
    # Make sure the bit representation is enumerated along a new axis, the
    #   very last axis
    axis = x.ndim
    # unpackbits() returns the bit representation in big-endian order, so we
    #   flip the array (with -8) to get litte-endian order
    return np.unpackbits(x[...,None], axis = axis)[...,-8:]

Unpacks an arbitrary decimal NumPy array into a binary representation along a new axis. Assumes decimal digits are on the interval [0, 255], i.e., that only 8-bit representations are needed.

Parameters

x : numpy.ndarray
 

Returns

numpy.ndarray
 
def geotransform(hdf,
ps=463.31271653,
nrow=2400,
ncol=2400,
metadata=re.compile('.*XDim=(?P<xdim>\\d+).*YDim=(?P<ydim>\\d+).*UpperLeftPointMtrs=\\((?P<ul>[0-9,\\-\\.]+)\\).*LowerRightMtrs=\\((?P<lr>[0-9,\\-\\.]+)\\).*', re.DOTALL))
Expand source code
def geotransform(
        hdf, ps = 463.31271653, nrow = 2400, ncol = 2400,
        metadata = VIIRS_METADATA):
    '''
    Prescribe a geotransform tuple for the output GeoTIFF. For MODIS/VIIRS
    sinsuoidal projections, the lower right corner coordinates are "the only
    metadata that accurately reflect the extreme corners of the gridded image"
    (Myneni et al. 2018, VIIRS LAI/fPAR User Guide). So, instead of using the
    reported upper-left (UL) corner coordinates, it is more accurate to use
    the lower-right (LR) corner coordinates and calculate the position of the
    UL corner based on the width and height of the image and the pixel size.
    NOTE that a rather odd pixel size is required to get the correct results
    verified against the HDF-EOS-to-GeoTIFF (HEG) Conversion Tool; see also
    Giglio et al. (2018), "Collection 6 MODIS Burned Area Product User's
    Guide, Version 1" Table 1.

        https://modis-land.gsfc.nasa.gov/pdf/MODIS_C6_BA_User_Guide_1.2.pdf

    Parameters
    ----------
    hdf : h5py.File
    ps : int or float
        The pixel size; in units matching the linear units of the SRS
        (Default: 463.3127 meters)
    nrow : int
        Number of rows in the output image (Default: 2400 for MODIS/VIIRS)
    ncol : int
        Number of columns in the output image (Default: 2400 for MODIS/VIIRS)
    metadata : re.Pattern
        Compiled regex that captures important metadata fields

    Returns
    -------
    tuple
        (x_min, pixel_width, 0, y_max, 0, -pixel_height)
    '''
    meta = hdf['HDFEOS INFORMATION/StructMetadata.0'][()].decode()
    lr = VIIRS_METADATA.match(meta).groupdict()['lr'].split(',')
    return ( # Subtract distance (meters) from LR corner to obtain UR corner
        float(lr[0]) - (ncol * ps), ps, 0, float(lr[1]) + (nrow * ps), 0, -ps)

Prescribe a geotransform tuple for the output GeoTIFF. For MODIS/VIIRS sinsuoidal projections, the lower right corner coordinates are "the only metadata that accurately reflect the extreme corners of the gridded image" (Myneni et al. 2018, VIIRS LAI/fPAR User Guide). So, instead of using the reported upper-left (UL) corner coordinates, it is more accurate to use the lower-right (LR) corner coordinates and calculate the position of the UL corner based on the width and height of the image and the pixel size. NOTE that a rather odd pixel size is required to get the correct results verified against the HDF-EOS-to-GeoTIFF (HEG) Conversion Tool; see also Giglio et al. (2018), "Collection 6 MODIS Burned Area Product User's Guide, Version 1" Table 1.

<https://modis-land.gsfc.nasa.gov/pdf/MODIS_C6_BA_User_Guide_1.2.pdf>

Parameters

hdf : h5py.File
 
ps : int or float
The pixel size; in units matching the linear units of the SRS (Default: 463.3127 meters)
nrow : int
Number of rows in the output image (Default: 2400 for MODIS/VIIRS)
ncol : int
Number of columns in the output image (Default: 2400 for MODIS/VIIRS)
metadata : re.Pattern
Compiled regex that captures important metadata fields

Returns

tuple
(x_min, pixel_width, 0, y_max, 0, -pixel_height)
def mod15a2h_qc_fail(x)
Expand source code
def mod15a2h_qc_fail(x):
    '''
    Returns pass/fail for QC flags based on the L4C fPAR QC protocol for the
    `FparLai_QC` band: Bad pixels have either `1` in the first bit ("Pixel not
    produced at all") or anything other than `00` ("clear") in bits 3-4.
    Output array is True wherever the array fails QC criteria. Compare to:

        np.vectorize(lambda v: v[0] == 1 or v[3:5] != '00')

    Parameters
    ----------
    x : numpy.ndarray
        Array where the last axis enumerates the unpacked bits
        (ones and zeros)

    Returns
    -------
    numpy.ndarray
        Boolean array with True wherever QC criteria are failed
    '''
    y = dec2bin_unpack(x)
    # Emit 1 = FAIL if these two bits are not == "00"
    c1 = y[...,3:5].sum(axis = -1).astype(np.uint8)
    # Emit 1 = FAIL if 1st bit == 1 ("Pixel not produced at all")
    c2 = y[...,0]
    # Intermediate arrays are 1 = FAIL, 0 = PASS
    return (c1 + c2) > 0

Returns pass/fail for QC flags based on the L4C fPAR QC protocol for the FparLai_QC band: Bad pixels have either 1 in the first bit ("Pixel not produced at all") or anything other than 00 ("clear") in bits 3-4. Output array is True wherever the array fails QC criteria. Compare to:

np.vectorize(lambda v: v[0] == 1 or v[3:5] != '00')

Parameters

x : numpy.ndarray
Array where the last axis enumerates the unpacked bits (ones and zeros)

Returns

numpy.ndarray
Boolean array with True wherever QC criteria are failed
def modis_from_wgs84(coords)
Expand source code
def modis_from_wgs84(coords):
    '''
    Given longitude-latitude coordinates, return the coordinates on the
    sinusoidal projection plane.

    Parameters
    ----------
    coords : tuple or list
        (Longitude, Latitude) coordinate pair

    Returns
    -------
    tuple
        (X, Y) coordinate pair in MODIS sinusoidal projection
    '''
    x, y = map(np.deg2rad, coords)
    return (SPHERE_RADIUS * x * np.cos(y), SPHERE_RADIUS * y)

Given longitude-latitude coordinates, return the coordinates on the sinusoidal projection plane.

Parameters

coords : tuple or list
(Longitude, Latitude) coordinate pair

Returns

tuple
(X, Y) coordinate pair in MODIS sinusoidal projection
def modis_row_col_from_wgs84(coords, nominal=500)
Expand source code
def modis_row_col_from_wgs84(coords, nominal = 500):
    '''
    Given longitude-latitude coordinates, return the corresponding row-column
    coordinates within a MODIS tile. NOTE: You'll need to determine which
    MODIS tile contains this row-column index with `modis_tile_from_wgs84()`.

    Parameters
    ----------
    coords : tuple or list
        (X, Y) coordinate pair in WGS84 (longitude, latitude)
    nominal : int
        Nominal resolution of MODIS raster: 250 (meters), 500, or 1000

    Returns
    -------
    tuple
        (Row, Column) coordinates
    '''
    x, y = modis_from_wgs84(coords) # Get coordinates in the projection plane
    # Get actual size of, e.g., "500-m" MODIS sinusoidal cell
    res = TILE_SIZE / float(TILE_LINES[nominal])
    return (
        np.floor((((YMAX - y) % TILE_SIZE) / res) - 0.5),
        np.floor((((x - XMIN) % TILE_SIZE) / res) - 0.5),
    )

Given longitude-latitude coordinates, return the corresponding row-column coordinates within a MODIS tile. NOTE: You'll need to determine which MODIS tile contains this row-column index with modis_tile_from_wgs84().

Parameters

coords : tuple or list
(X, Y) coordinate pair in WGS84 (longitude, latitude)
nominal : int
Nominal resolution of MODIS raster: 250 (meters), 500, or 1000

Returns

tuple
(Row, Column) coordinates
def modis_row_col_to_wgs84(coords, h, v, nominal=500)
Expand source code
def modis_row_col_to_wgs84(coords, h, v, nominal = 500):
    '''
    Convert pixel coordinates in a specific MODIS tile to longitude-latitude
    coordinates.

    Parameters
    ----------
    coords : tuple or list
        (X, Y) coordinate pair in MODIS sinusoidal projection
    h : int
        MODIS tile "h" index
    v : int
        MODIS tile "v" index
    nominal : int
        Nominal resolution of MODIS raster: 250 (meters), 500, or 1000

    Returns
    -------
    tuple
        (Longitude, Latitude) coordinates
    '''
    r, c = coords
    lines = TILE_LINES[nominal]
    assert (0 <= r <= (lines - 1)) and (0 <= c <= (lines - 1)),\
        'Row and column indices must be in the range [0, %d]' % (lines - 1)
    # Get actual size of, e.g., "500-m" MODIS sinusoidal cell
    res = TILE_SIZE / float(TILE_LINES[nominal])
    x = ((c + 0.5) * res) + (h * TILE_SIZE) + XMIN
    y = YMAX - ((r + 0.5) * res) - (v * TILE_SIZE)
    return modis_to_wgs84((x, y))

Convert pixel coordinates in a specific MODIS tile to longitude-latitude coordinates.

Parameters

coords : tuple or list
(X, Y) coordinate pair in MODIS sinusoidal projection
h : int
MODIS tile "h" index
v : int
MODIS tile "v" index
nominal : int
Nominal resolution of MODIS raster: 250 (meters), 500, or 1000

Returns

tuple
(Longitude, Latitude) coordinates
def modis_tile_from_wgs84(coords)
Expand source code
def modis_tile_from_wgs84(coords):
    '''
    Given longitude-latitude coordinates, return the MODIS tile (H,V) that
    contains them.

    Parameters
    ----------
    coords : tuple or list
        (Longitude, Latitude) coordinate pair

    Returns
    -------
    tuple
        (H,V) tile identifier
    '''
    x, y = modis_from_wgs84(coords) # Get coordinates in the projection plane
    return (
        np.floor((x - XMIN) / TILE_SIZE),
        np.floor((YMAX - y) / TILE_SIZE))

Given longitude-latitude coordinates, return the MODIS tile (H,V) that contains them.

Parameters

coords : tuple or list
(Longitude, Latitude) coordinate pair

Returns

tuple
(H,V) tile identifier
def modis_to_wgs84(coords)
Expand source code
def modis_to_wgs84(coords):
    '''
    Convert coordinates on the MODIS sinusoidal plane to longitude-latitude
    coordinates (WGS84).

    Parameters
    ----------
    coords : tuple or list
        (X, Y) coordinate pair in MODIS sinusoidal projection

    Returns
    -------
    tuple
        (Longitude, Latitude) coordinate pair
    '''
    x, y = coords
    lat = y / SPHERE_RADIUS
    lng = x / (SPHERE_RADIUS * np.cos(lat))
    return tuple(map(np.rad2deg, (lng, lat)))

Convert coordinates on the MODIS sinusoidal plane to longitude-latitude coordinates (WGS84).

Parameters

coords : tuple or list
(X, Y) coordinate pair in MODIS sinusoidal projection

Returns

tuple
(Longitude, Latitude) coordinate pair
def vnp15a2h_cloud_fail(x)
Expand source code
def vnp15a2h_cloud_fail(x):
    '''
    Returns pass/fail for QC flags based on the L4C fPAR QC protocol for the
    `FparExtra_QC` band (cloud QC band): Bad pixels have anything OTHER THAN
    `1` second least-significant bit; `00` and `01` being acceptable cloud QC
    flags ("Confident clear" and "Probably clear", respectively).

    Parameters
    ----------
    x : numpy.ndarray
        Unsigned, 8-bit integer array

    Returns
    -------
    numpy.ndarray
        Boolean array
    '''
    y = dec2bin_unpack(x)
    return y[...,-2] > 0

Returns pass/fail for QC flags based on the L4C fPAR QC protocol for the FparExtra_QC band (cloud QC band): Bad pixels have anything OTHER THAN 1 second least-significant bit; 00 and 01 being acceptable cloud QC flags ("Confident clear" and "Probably clear", respectively).

Parameters

x : numpy.ndarray
Unsigned, 8-bit integer array

Returns

numpy.ndarray
Boolean array
def vnp15a2h_qc_fail(x)
Expand source code
def vnp15a2h_qc_fail(x):
    '''
    Returns pass/fail for QC flags based on the L4C fPAR QC protocol for the
    `FparLai_QC` band: Bad pixels have either `11` in the first two bits
    ("Fill Value") or anything other than `0` in the 3rd least-significant
    bits, which combines "Pixel not produced at all". For example, see decimal
    number 80:

        0101|0|000 where "000" is the combined (Fill bit | Retrieval quality)

    Parameters
    ----------
    x : numpy.ndarray
        Unsigned, 8-bit integer array

    Returns
    -------
    numpy.ndarray
        Boolean array
    '''
    y = dec2bin_unpack(x)
    # Emit 1 = FAIL if sum("11") == 2; "BiomeType" == "Filled Value"
    c1 = np.where(y[...,0:2].sum(axis = -1) == 2, 1, 0).astype(np.uint8)
    # Emit 1 = FAIL if 3rd bit == 1; "SCF_QC" == "Pixel not produced at all"
    c2 = y[...,5]
    # Intermediate arrays are 1 = FAIL, 0 = PASS
    return (c1 + c2) > 0

Returns pass/fail for QC flags based on the L4C fPAR QC protocol for the FparLai_QC band: Bad pixels have either 11 in the first two bits ("Fill Value") or anything other than 0 in the 3rd least-significant bits, which combines "Pixel not produced at all". For example, see decimal number 80:

0101|0|000 where "000" is the combined (Fill bit | Retrieval quality)

Parameters

x : numpy.ndarray
Unsigned, 8-bit integer array

Returns

numpy.ndarray
Boolean array