Analysis: trends and summaries over raster stacks on mixed grids and over multi-band GeoTIFFs dated by their band descriptions (#579)

crates/data/src/raster/geotiff.rs

@@ -878,6 +878,13 @@ pub(super) fn read_limited(r: Reader<'_>, name: &str, max: usize, open: Option<O
     let mut dec = Decoder::new(r).map_err(err)?.with_limits(limits.clone());
     let (full_w, full_h) = dec.dimensions().map_err(err)?;
     let geo = geo_tags(&mut dec)?;
+    // The bands' descriptions GDAL keeps in the file (an `.aux.xml` sidecar's are applied after).
+    let described = dec
+        .find_tag(Tag::Unknown(42112))
+        .ok()
+        .flatten()
+        .and_then(|v| v.into_string().ok())
+        .and_then(|t| super::pam::parse_gdal_metadata(&t));
     let spp = samples(&mut dec);
     let bytes = sample_bytes(&mut dec);
     // The decoded image is held at its own pixel size before conversion to f32, so bound bytes
@@ -1041,6 +1048,9 @@ pub(super) fn read_limited(r: Reader<'_>, name: &str, max: usize, open: Option<O
         let mut r = Raster::new(name, w, h, bands, super::unreferenced(w, h), Crs::LonLat, geo.nodata, color, data_type)?;
         r.georeferenced = false;
         r.overview = overview;
+        if let Some(p) = &described {
+            r.apply_pam(p);
+        }
         return Ok(r);
     }
     let mut transform = geo.transform;
@@ -1087,6 +1097,9 @@ pub(super) fn read_limited(r: Reader<'_>, name: &str, max: usize, open: Option<O
     };
     let mut r = Raster::new(name, w, h, bands, transform, crs, geo.nodata, color, data_type)?;
     r.overview = overview;
+    if let Some(p) = &described {
+        r.apply_pam(p);
+    }
     Ok(r)
 }
 
@@ -1906,6 +1919,29 @@ mod tests {
         assert!(read.as_ref().err().is_none_or(|e| !e.to_string().contains("fewer than 8 bits")), "{read:?}");
     }
 
+    #[test]
+    fn band_descriptions_in_the_file_name_its_bands() {
+        // GDAL and rasterio keep band descriptions in the GDAL_METADATA tag (42112).
+        use tiff::encoder::{TiffEncoder, colortype};
+        let mut out = std::io::Cursor::new(Vec::new());
+        {
+            let mut enc = TiffEncoder::new(&mut out).unwrap();
+            let mut img = enc.new_image::<colortype::Gray32Float>(2, 2).unwrap();
+            let e = img.encoder();
+            e.write_tag(Tag::ModelPixelScaleTag, &[0.5f64, 0.5, 0.0][..]).unwrap();
+            e.write_tag(Tag::ModelTiepointTag, &[0.0f64, 0.0, 0.0, 10.0, 50.0, 0.0][..]).unwrap();
+            e.write_tag(Tag::GeoKeyDirectoryTag, &[1u16, 1, 0, 2, 1024, 0, 1, 2, 2048, 0, 1, 4326][..]).unwrap();
+            let xml = r#"<GDALMetadata><Item name="DESCRIPTION" sample="0" role="description">NDVI 2024-06-01</Item></GDALMetadata>"#;
+            e.write_tag(Tag::Unknown(42112), xml).unwrap();
+            img.write_data(&[0.1f32, 0.2, 0.3, 0.4]).unwrap();
+        }
+        let r = read(std::io::Cursor::new(out.into_inner()), "ndvi").unwrap();
+        assert_eq!(r.bands[0].name, "NDVI 2024-06-01");
+        // Without the tag, the band keeps its number.
+        let plain = read(std::io::Cursor::new(grey8(2, 2, &[1, 2, 3, 4], None)), "plain").unwrap();
+        assert_ne!(plain.bands[0].name, "NDVI 2024-06-01");
+    }
+
     #[test]
     fn a_files_own_projection_is_read_and_written_back() {
         // The NLCD's Albers, as its GeoKeys give it.

crates/data/src/raster/pam.rs

@@ -118,6 +118,42 @@ pub fn parse(text: &str) -> Option<Pam> {
     Some(pam)
 }
 
+/// The band descriptions a GeoTIFF holds itself, in GDAL's metadata tag (`GDAL_METADATA`, 42112),
+/// where GDAL and rasterio write them: `None` when it isn't such a document or describes no band.
+/// Its `sample` counts from 0.
+///
+/// ```xml
+/// <GDALMetadata>
+///   <Item name="DESCRIPTION" sample="0" role="description">NDVI 2024-06-01</Item>
+/// </GDALMetadata>
+/// ```
+pub fn parse_gdal_metadata(text: &str) -> Option<Pam> {
+    let opts = ParsingOptions { allow_dtd: false, nodes_limit: crate::io::XML_NODES, ..ParsingOptions::default() };
+    let doc = crate::xml::parse(text, opts).ok()?;
+    let root = doc.root_element();
+    if root.tag_name().name() != "GDALMetadata" {
+        return None;
+    }
+    let mut pam = Pam::default();
+    // The bands described so far, so a long document isn't searched again for each item.
+    let mut seen = std::collections::HashSet::new();
+    for item in root.children().filter(|n| n.has_tag_name("Item") && n.attribute("role") == Some("description")) {
+        let band = item
+            .attribute("sample")
+            .and_then(|v| v.trim().parse::<usize>().ok())
+            .and_then(|s| s.checked_add(1))
+            .filter(|n| *n <= MAX_BANDS);
+        let description = item.text().map(str::trim).filter(|d| !d.is_empty()).map(|d| d.chars().take(256).collect::<String>());
+        // A band's first description, if it is given twice.
+        if let (Some(band), Some(description)) = (band, description)
+            && seen.insert(band)
+        {
+            pam.bands.push(PamBand { band, description: Some(description), ..PamBand::default() });
+        }
+    }
+    (!pam.bands.is_empty()).then_some(pam)
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -164,6 +200,29 @@ mod tests {
         assert!(parse(&deep).is_none());
     }
 
+    #[test]
+    fn reads_band_descriptions_from_gdal_metadata() {
+        let xml = r#"<GDALMetadata>
+  <Item name="AREA_OR_POINT">Area</Item>
+  <Item name="DESCRIPTION" sample="1" role="description"> NDVI 2024-07-01 </Item>
+  <Item name="DESCRIPTION" sample="0" role="description">NDVI 2024-06-01</Item>
+  <Item name="DESCRIPTION" sample="0" role="description">again</Item>
+  <Item name="SCALE" sample="0" role="scale">0.0001</Item>
+  <Item name="DESCRIPTION" sample="x" role="description">no band</Item>
+  <Item name="DESCRIPTION" sample="18446744073709551615" role="description">past the last band</Item>
+  <Item name="DESCRIPTION" sample="65536" role="description">past the last band</Item>
+  <Item name="DESCRIPTION" sample="2" role="description"></Item>
+</GDALMetadata>"#;
+        let pam = parse_gdal_metadata(xml).unwrap();
+        assert_eq!(pam.bands.len(), 2, "{pam:?}");
+        assert_eq!(pam.band(1).unwrap().description.as_deref(), Some("NDVI 2024-06-01"));
+        assert_eq!(pam.band(2).unwrap().description.as_deref(), Some("NDVI 2024-07-01"));
+        assert_eq!(pam.nodata(), None);
+        for text in ["", "not xml", "<GDALMetadata/>", "<PAMDataset/>", "<GDALMetadata><Item name=\"A\">1</Item></GDALMetadata>"] {
+            assert!(parse_gdal_metadata(text).is_none(), "{text}");
+        }
+    }
+
     #[test]
     fn sidecar_names() {
         let p = sidecar_paths(Path::new("/d/scene.tif"));