From d86f9787383d9384fac6cb0a96958a3f894c605a Mon Sep 17 00:00:00 2001 From: Chuck Daniels Date: Fri, 14 Jun 2024 19:05:15 -0400 Subject: [PATCH] Correct VI metadata Fixes #11 --- hls_vi/generate_indices.py | 108 +++--- ...-VI.S30.T13RCN.2024128T173909.v2.0.cmr.xml | 321 ++++++++++++++++++ tests/test_vi.py | 37 +- 3 files changed, 410 insertions(+), 56 deletions(-) create mode 100644 tests/fixtures/HLS-VI.S30.T13RCN.2024128T173909.v2.0.cmr.xml diff --git a/hls_vi/generate_indices.py b/hls_vi/generate_indices.py index 4c1cf74..2addffe 100644 --- a/hls_vi/generate_indices.py +++ b/hls_vi/generate_indices.py @@ -23,6 +23,26 @@ BandData: TypeAlias = Mapping["Band", np.ma.masked_array] IndexFunction = Callable[[BandData], np.ma.masked_array] +fixed_tags = ( + "add_offset", + "ACCODE", + "AREA_OR_POINT", + "cloud_coverage", + "HORIZONTAL_CS_NAME", + "MEAN_SUN_AZIMUTH_ANGLE", + "MEAN_SUN_ZENITH_ANGLE", + "MEAN_VIEW_AZIMUTH_ANGLE", + "MEAN_VIEW_ZENITH_ANGLE", + "NBAR_SOLAR_ZENITH", + "NCOLS", + "NROWS", + "SENSING_TIME", + "spatial_coverage", + "SPATIAL_RESOLUTION", + "ULX", + "ULY", +) + @unique class Band(Enum): @@ -59,11 +79,20 @@ class S30Band(InstrumentBand): class Instrument(Enum): - L30 = L30Band - S30 = S30Band - - def __init__(self, band_type: Type[InstrumentBand]) -> None: + # Example LANDSAT_PRODUCT_ID: LC08_L1TP_069014_20240429_20240430_02_RT + # Pull out satellite number ("08"), convert to an integer, and prepend "L" to get + # the satellite name. + L30 = L30Band, lambda tags: f"L{int(tags.get('LANDSAT_PRODUCT_ID', '')[2:4])}" + # Example PRODUCT_URI: S2B_MSIL1C_... + # Simply pull off the first part of the string before the first underscore to get + # the satellite name. + S30 = S30Band, lambda tags: tags.get("PRODUCT_URI", "").split("_")[0] + + def __init__( + self, band_type: Type[InstrumentBand], parse_satellite: Callable[[Tags], str] + ) -> None: self.bands = list(band_type) + self.parse_satellite = parse_satellite @classmethod def named(cls, name: str) -> "Instrument": @@ -73,6 +102,9 @@ def named(cls, name: str) -> "Instrument": raise ValueError(f"Invalid instrument name: {name}") + def satellite(self, tags: Tags) -> str: + return self.parse_satellite(tags) + @dataclass class GranuleId: @@ -137,7 +169,7 @@ def read_granule_bands(input_dir: Path) -> Granule: with rasterio.open(input_dir / tifnames[0]) as tif: crs = tif.crs transform = tif.transform - tags = select_tags(tif.tags()) + tags = select_tags(id_, tif.tags()) return Granule(id_, crs, transform, tags, dict(zip(harmonized_bands, data))) @@ -158,7 +190,7 @@ def apply_fmask(data: np.ndarray, fmask: np.ndarray) -> np.ma.masked_array: return np.ma.masked_array(data, fmask & cloud_like != 0) -def select_tags(tags: Tags) -> Tags: +def select_tags(granule_id: GranuleId, tags: Tags) -> Tags: """ Selects tags from the input tags that are relevant to the HLS VI product. @@ -169,19 +201,9 @@ def select_tags(tags: Tags) -> Tags: Mapping of relevant VI tags. """ return { - "ACCODE": tags["ACCODE"], - "cloud_coverage": tags.get("cloud_coverage"), - "HORIZONTAL_CS_NAME": tags.get("HORIZONTAL_CS_NAME"), - "MEAN_SUN_AZIMUTH_ANGLE": tags.get("MEAN_SUN_AZIMUTH_ANGLE"), - "MEAN_SUN_ZENITH_ANGLE": tags.get("MEAN_SUN_ZENITH_ANGLE"), - "MEAN_VIEW_AZIMUTH_ANGLE": tags.get("MEAN_VIEW_AZIMUTH_ANGLE"), - "MEAN_VIEW_ZENITH_ANGLE": tags.get("MEAN_VIEW_ZENITH_ANGLE"), - "NBAR_SOLAR_ZENITH": tags.get("NBAR_SOLAR_ZENITH"), - "SPACECRAFT_NAME": tags.get("SPACECRAFT_NAME"), - "TILE_ID": tags.get("SENTINEL2_TILEID"), - "SENSING_TIME": tags.get("SENSING_TIME"), - "SENSOR": tags.get("SENSOR"), - "spatial_coverage": tags.get("spatial_coverage"), + **{tag: tags.get(tag) for tag in fixed_tags}, + "MGRS_TILE_ID": granule_id.tile_id, + "SATELLITE": granule_id.instrument.satellite(tags), } @@ -231,8 +253,10 @@ def write_granule_index( dst.write(data.filled(), 1) dst.update_tags( **granule.tags, - longname=index.value, + long_name=index.long_name, + scale_factor=index.scale_factor, HLS_VI_PROCESSING_TIME=processing_time, + _FillValue=data.fill_value, ) # Create browse image using NDVI @@ -242,14 +266,14 @@ def write_granule_index( def evi(data: BandData) -> np.ma.masked_array: b, r, nir = data[Band.B], data[Band.R], data[Band.NIR] - return 10_000 * 2.5 * (nir - r) / (nir + 6 * r - 7.5 * b + 1) # type: ignore + return 2.5 * (nir - r) / (nir + 6 * r - 7.5 * b + 1) # type: ignore def msavi(data: BandData) -> np.ma.masked_array: r, nir = data[Band.R], data[Band.NIR] sqrt_term = (2 * nir + 1) ** 2 - 8 * (nir - r) # type: ignore - result: np.ma.masked_array = 10_000 * np.ma.where( + result: np.ma.masked_array = np.ma.where( sqrt_term >= 0, (2 * nir + 1 - np.sqrt(sqrt_term)) / 2, # type: ignore np.nan, @@ -261,32 +285,32 @@ def msavi(data: BandData) -> np.ma.masked_array: def nbr(data: BandData) -> np.ma.masked_array: nir, swir2 = data[Band.NIR], data[Band.SWIR2] - return 10_000 * (nir - swir2) / (nir + swir2) # type: ignore + return (nir - swir2) / (nir + swir2) # type: ignore def nbr2(data: BandData) -> np.ma.masked_array: swir1, swir2 = data[Band.SWIR1], data[Band.SWIR2] - return 10_000 * (swir1 - swir2) / (swir1 + swir2) # type: ignore + return (swir1 - swir2) / (swir1 + swir2) # type: ignore def ndmi(data: BandData) -> np.ma.masked_array: nir, swir1 = data[Band.NIR], data[Band.SWIR1] - return 10_000 * (nir - swir1) / (nir + swir1) # type: ignore + return (nir - swir1) / (nir + swir1) # type: ignore def ndvi(data: BandData) -> np.ma.masked_array: r, nir = data[Band.R], data[Band.NIR] - return 10_000 * (nir - r) / (nir + r) # type: ignore + return (nir - r) / (nir + r) # type: ignore def ndwi(data: BandData) -> np.ma.masked_array: g, nir = data[Band.G], data[Band.NIR] - return 10_000 * (g - nir) / (g + nir) # type: ignore + return (g - nir) / (g + nir) # type: ignore def savi(data: BandData) -> np.ma.masked_array: r, nir = data[Band.R], data[Band.NIR] - return 10_000 * 1.5 * (nir - r) / (nir + r + 0.5) # type: ignore + return 1.5 * (nir - r) / (nir + r + 0.5) # type: ignore def tvi(data: BandData) -> np.ma.masked_array: @@ -296,28 +320,30 @@ def tvi(data: BandData) -> np.ma.masked_array: class Index(Enum): - EVI = "Enhanced Vegetation Index" - MSAVI = "Modified Soil-Adjusted Vegetation Index" - NBR = "Normalized Burn Ratio" - NBR2 = "Normalized Burn Ratio 2" - NDMI = "Normalized Difference Moisture Index" - NDVI = "Normalized Difference Vegetation Index" - NDWI = "Normalized Difference Water Index" - SAVI = "Soil-Adjusted Vegetation Index" - TVI = "Triangular Vegetation Index" - - def __init__(self, longname: str) -> None: + EVI = ("Enhanced Vegetation Index",) + MSAVI = ("Modified Soil-Adjusted Vegetation Index",) + NBR = ("Normalized Burn Ratio",) + NBR2 = ("Normalized Burn Ratio 2",) + NDMI = ("Normalized Difference Moisture Index",) + NDVI = ("Normalized Difference Vegetation Index",) + NDWI = ("Normalized Difference Water Index",) + SAVI = ("Soil-Adjusted Vegetation Index",) + TVI = ("Triangular Vegetation Index", 1) + + def __init__(self, long_name: str, scale_factor: float = 0.0001) -> None: function_name = self.name.lower() index_function: Optional[IndexFunction] = globals().get(function_name) if not index_function or not callable(index_function): raise ValueError(f"Index function not found: {function_name}") - self.longname = longname + self.long_name = long_name self.compute_index = index_function + self.scale_factor = scale_factor def __call__(self, data: BandData) -> np.ma.masked_array: - return np.ma.round(self.compute_index(data)).astype(np.int16) + scaled_index = self.compute_index(data) / self.scale_factor + return np.ma.round(scaled_index).astype(np.int16) def parse_args() -> Tuple[Path, Path]: diff --git a/tests/fixtures/HLS-VI.S30.T13RCN.2024128T173909.v2.0.cmr.xml b/tests/fixtures/HLS-VI.S30.T13RCN.2024128T173909.v2.0.cmr.xml new file mode 100644 index 0000000..45d76fa --- /dev/null +++ b/tests/fixtures/HLS-VI.S30.T13RCN.2024128T173909.v2.0.cmr.xml @@ -0,0 +1,321 @@ + + + HLS-VI.S30.T13RCN.2024128T173909.v2.0 + 2024-05-09T20:08:04.550553Z + 2024-05-09T20:08:04.550572Z + + Update HLS-VI New Collection ID String + + + HLS-VI.S30.T13RCN.2024128T173909 + DAY + UPDATE HLS Prodution DATETIME + 2.0 + + + + 2024-05-09T15:07:56.146630Z + 2024-05-09T15:07:56.146630Z + + + + + + + + + -107.06959008 + 29.81425505 + + + -105.93358748 + 29.82716844 + + + -105.92460396 + 28.8362971 + + + -107.04968888 + 28.82389371 + + + + + + + + + Sentinel-2B + + + Sentinel-2 MSI + + + + + + + PRODUCT_URI + + S2B_MSIL1C_20240507T173909_N0510_R098_T13RCN_20240507T205106.SAFE + + + + CLOUD_COVERAGE + + 9 + + + + MGRS_TILE_ID + + 13RCN + + + + SPATIAL_COVERAGE + + 100 + + + + SPATIAL_RESOLUTION + + 30.0 + + + + HLS_PROCESSING_TIME + + 2024-05-09T13:06:35Z + + + + SENSING_TIME + + 2024-05-07T17:55:57.208242Z + + + + HORIZONTAL_CS_CODE + + EPSG:32613 + + + + HORIZONTAL_CS_NAME + + WGS84 / UTM zone 13N + + + + ULX + + 300000.0 + + + + ULY + + 3300000.0 + + + + SPATIAL_RESAMPLING_ALG + + Area Weighted Average + + + + ADD_OFFSET + + 0 + + + + REF_SCALE_FACTOR + + 0.0001 + + + + ANG_SCALE_FACTOR + + 0.01 + + + + FILLVALUE + + -9999 + + + + QA_FILLVALUE + + 255 + + + + NCOLS + + 3660 + + + + NROWS + + 3660 + + + + MEAN_SUN_AZIMUTH_ANGLE + + 125.08344994 + + + + MEAN_SUN_ZENITH_ANGLE + + 19.51903921 + + + + MEAN_VIEW_AZIMUTH_ANGLE + + 198.11772041 + + + + MEAN_VIEW_ZENITH_ANGLE + + 2.85707672 + + + + NBAR_SOLAR_ZENITH + + 21.48626569 + + + + MSI_BAND_01_BANDPASS_ADJUSTMENT_SLOPE_AND_OFFSET + + 0.9959 + -0.0002 + + + + MSI_BAND_02_BANDPASS_ADJUSTMENT_SLOPE_AND_OFFSET + + 0.9778 + -0.004 + + + + MSI_BAND_03_BANDPASS_ADJUSTMENT_SLOPE_AND_OFFSET + + 1.0075 + -0.0008 + + + + MSI_BAND_04_BANDPASS_ADJUSTMENT_SLOPE_AND_OFFSET + + 0.9761 + 0.001 + + + + MSI_BAND_11_BANDPASS_ADJUSTMENT_SLOPE_AND_OFFSET + + 1.0 + -0.0003 + + + + MSI_BAND_12_BANDPASS_ADJUSTMENT_SLOPE_AND_OFFSET + + 0.9867 + 0.0004 + + + + MSI_BAND_8A_BANDPASS_ADJUSTMENT_SLOPE_AND_OFFSET + + 0.9966 + 0.0 + + + + AROP_AVE_XSHIFT(METERS) + + 0.0 + + + + AROP_AVE_YSHIFT(METERS) + + 0.0 + + + + AROP_NCP + + 0 + + + + AROP_RMSE(METERS) + + 0.0 + + + + AROP_S2_REFIMG + + NONE + + + + ACCODE + + LaSRC v3.2.0 + + + + PROCESSING_BASELINE + + 05.10 + + + + IDENTIFIER_PRODUCT_DOI + + 10.5067/HLS/HLSS30.002 + + + + IDENTIFIER_PRODUCT_DOI_AUTHORITY + + https://doi.org + + + + + + + + Cloud Optimized GeoTIFF (COG) + + + https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-public/HLSS30.020/HLS.S30.T13RCN.2024128T173909.v2.0/HLS.S30.T13RCN.2024128T173909.v2.0.jpg + Download HLS.S30.T13RCN.2024128T173909.v2.0.jpg + + + s3://lp-prod-public/HLSS30.020/HLS.S30.T13RCN.2024128T173909.v2.0/HLS.S30.T13RCN.2024128T173909.v2.0.jpg + This link provides direct download access via S3 to the granule + + + diff --git a/tests/test_vi.py b/tests/test_vi.py index cdb8b55..812ac4b 100644 --- a/tests/test_vi.py +++ b/tests/test_vi.py @@ -1,5 +1,6 @@ from pathlib import Path from xml.etree import ElementTree as ET +import contextlib import io import os @@ -26,18 +27,25 @@ def tifs_equal(tif1: Path, tif2: Path): ) -def remove_datetime_elements(tree: ET.ElementTree) -> ET.ElementTree: - root = tree.getroot() +def remove_element(root: ET.Element, path: str) -> None: + parent_path = "/".join(path.split("/")[:-1]) + parent = root.find(parent_path) + child = root.find(path) + + assert parent is not None + assert child is not None + + parent.remove(child) - root.remove(tree.find("./InsertTime")) - root.remove(tree.find("./LastUpdate")) - data_granule = root.find("./DataGranule") - data_granule.remove(data_granule.find("./ProductionDateTime")) +def remove_datetime_elements(tree: ET.ElementTree) -> ET.ElementTree: + root = tree.getroot() - range_date_time = root.find("./Temporal/RangeDateTime") - range_date_time.remove(range_date_time.find("./BeginningDateTime")) - range_date_time.remove(range_date_time.find("./EndingDateTime")) + remove_element(root, "./InsertTime") + remove_element(root, "./LastUpdate") + remove_element(root, "./DataGranule/ProductionDateTime") + remove_element(root, "./Temporal/RangeDateTime/BeginningDateTime") + remove_element(root, "./Temporal/RangeDateTime/EndingDateTime") return tree @@ -73,11 +81,10 @@ def test_generate_indices(input_dir, tmp_path: Path): "tests/fixtures/HLS.L30.T06WVS.2024120T211159.v2.0", "tests/fixtures/HLS-VI.L30.T06WVS.2024120T211159.v2.0", ), - # TODO: Create tests/fixtures/HLS-VI.S30.T13RCN.2024128T173909.v2.0/HLS-VI.S30.T13RCN.2024128T173909.v2.0.cmr.xml # noqa: E501 - # ( - # "tests/fixtures/HLS.S30.T13RCN.2024128T173909.v2.0", - # "tests/fixtures/HLS-VI.S30.T13RCN.2024128T173909.v2.0", - # ), + ( + "tests/fixtures/HLS.S30.T13RCN.2024128T173909.v2.0", + "tests/fixtures/HLS-VI.S30.T13RCN.2024128T173909.v2.0", + ), ], ) def test_generate_cmr_metadata(input_dir, output_dir): @@ -108,5 +115,5 @@ def test_generate_cmr_metadata(input_dir, output_dir): actual_metadata.getvalue().decode() == expected_metadata.getvalue().decode() ) finally: - if actual_metadata_path.exists(): + with contextlib.suppress(FileNotFoundError): actual_metadata_path.unlink()