Preserve GDALMetadata XML (tag 42112) through read/write

Band names, statistics, and other GDAL-specific metadata stored in
the GDALMetadata XML tag (42112) are now read, exposed, and written
back.

Read: the XML is parsed into a dict stored in attrs['gdal_metadata'].
Dataset-level items use string keys ('DataType'), per-band items use
(name, band_int) tuple keys (('STATISTICS_MAXIMUM', 0)). The raw XML
is also available in attrs['gdal_metadata_xml'].

Write: accepts gdal_metadata_xml on write(), or extracts from
DataArray attrs on write_geotiff(). If attrs has a gdal_metadata
dict but no raw XML, it's re-serialized automatically.

Round-trip verified on the USGS 1-meter DEM which has statistics,
layer type, and data type metadata -- all items survive intact.

7 new tests: XML parsing, dict serialization, file round-trip,
DataArray attrs preservation, no-metadata baseline, real-file read,
and real-file round-trip.

Author: brendancol

xrspatial/geotiff/__init__.py

@@ -176,6 +176,12 @@ def read_geotiff(source: str, *, window=None,
     if geo_info.vertical_units is not None:
         attrs['vertical_units'] = geo_info.vertical_units
 
+    # GDAL metadata (tag 42112)
+    if geo_info.gdal_metadata is not None:
+        attrs['gdal_metadata'] = geo_info.gdal_metadata
+    if geo_info.gdal_metadata_xml is not None:
+        attrs['gdal_metadata_xml'] = geo_info.gdal_metadata_xml
+
     # Resolution / DPI metadata
     if geo_info.x_resolution is not None:
         attrs['x_resolution'] = geo_info.x_resolution
@@ -275,6 +281,7 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
     x_res = None
     y_res = None
     res_unit = None
+    gdal_meta_xml = None
 
     # Resolve crs argument: can be int (EPSG) or str (WKT/PROJ)
     if isinstance(crs, int):
@@ -297,6 +304,13 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
             nodata = data.attrs.get('nodata')
         if data.attrs.get('raster_type') == 'point':
             raster_type = RASTER_PIXEL_IS_POINT
+        # GDAL metadata from attrs (prefer raw XML, fall back to dict)
+        gdal_meta_xml = data.attrs.get('gdal_metadata_xml')
+        if gdal_meta_xml is None:
+            gdal_meta_dict = data.attrs.get('gdal_metadata')
+            if isinstance(gdal_meta_dict, dict):
+                from ._geotags import _build_gdal_metadata_xml
+                gdal_meta_xml = _build_gdal_metadata_xml(gdal_meta_dict)
         # Resolution / DPI from attrs
         x_res = data.attrs.get('x_resolution')
         y_res = data.attrs.get('y_resolution')
@@ -326,6 +340,7 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
         x_resolution=x_res,
         y_resolution=y_res,
         resolution_unit=res_unit,
+        gdal_metadata_xml=gdal_meta_xml,
     )
 
 

xrspatial/geotiff/_geotags.py

@@ -109,10 +109,55 @@ class GeoInfo:
     vertical_units_code: int | None = None
     # WKT CRS string (resolved from EPSG via pyproj, or provided by caller)
     crs_wkt: str | None = None
+    # GDAL metadata: dict of {name: value} for dataset-level items,
+    # and {(name, band): value} for per-band items.  Raw XML also kept.
+    gdal_metadata: dict | None = None
+    gdal_metadata_xml: str | None = None
     # Raw geokeys dict for anything else
     geokeys: dict[int, int | float | str] = field(default_factory=dict)
 
 
+def _parse_gdal_metadata(xml_str: str) -> dict:
+    """Parse GDALMetadata XML into a flat dict.
+
+    Dataset-level items are stored as ``{name: value}``.
+    Per-band items are stored as ``{(name, band_int): value}``.
+    """
+    import xml.etree.ElementTree as ET
+    result = {}
+    try:
+        root = ET.fromstring(xml_str)
+        for item in root.findall('Item'):
+            name = item.get('name', '')
+            sample = item.get('sample')
+            text = item.text or ''
+            if sample is not None:
+                result[(name, int(sample))] = text
+            else:
+                result[name] = text
+    except ET.ParseError:
+        pass
+    return result
+
+
+def _build_gdal_metadata_xml(meta: dict) -> str:
+    """Serialize a metadata dict back to GDALMetadata XML.
+
+    Accepts the same dict format that _parse_gdal_metadata produces:
+    string keys for dataset-level, (name, band) tuples for per-band.
+    """
+    lines = ['<GDALMetadata>']
+    for key, value in meta.items():
+        if isinstance(key, tuple):
+            name, sample = key
+            lines.append(
+                f'  <Item name="{name}" sample="{sample}">{value}</Item>')
+        else:
+            lines.append(f'  <Item name="{key}">{value}</Item>')
+    lines.append('</GDALMetadata>')
+    return '\n'.join(lines) + '\n'
+
+
 def _epsg_to_wkt(epsg: int) -> str | None:
     """Resolve an EPSG code to a WKT string using pyproj.
 
@@ -379,6 +424,12 @@ def extract_geo_info(ifd: IFD, data: bytes | memoryview,
         except (ValueError, TypeError):
             pass
 
+    # Parse GDALMetadata XML (tag 42112)
+    gdal_metadata = None
+    gdal_metadata_xml = ifd.gdal_metadata
+    if gdal_metadata_xml is not None:
+        gdal_metadata = _parse_gdal_metadata(gdal_metadata_xml)
+
     # Extract palette colormap (Photometric=3, tag 320)
     colormap = None
     if ifd.photometric == 3:
@@ -430,6 +481,8 @@ def extract_geo_info(ifd: IFD, data: bytes | memoryview,
         vertical_units=vert_units_name,
         vertical_units_code=vert_units_code,
         crs_wkt=crs_wkt,
+        gdal_metadata=gdal_metadata,
+        gdal_metadata_xml=gdal_metadata_xml,
         geokeys=geokeys,
     )
 

xrspatial/geotiff/_header.py

@@ -35,6 +35,7 @@
 TAG_TILE_BYTE_COUNTS = 325
 TAG_COLORMAP = 320
 TAG_SAMPLE_FORMAT = 339
+TAG_GDAL_METADATA = 42112
 TAG_GDAL_NODATA = 42113
 
 # GeoTIFF tags
@@ -184,6 +185,16 @@ def colormap(self) -> tuple | None:
         """ColorMap tag (320) values, or None if absent."""
         return self.get_values(TAG_COLORMAP)
 
+    @property
+    def gdal_metadata(self) -> str | None:
+        """GDALMetadata XML string (tag 42112), or None if absent."""
+        v = self.get_value(TAG_GDAL_METADATA)
+        if v is None:
+            return None
+        if isinstance(v, bytes):
+            return v.rstrip(b'\x00').decode('ascii', errors='replace')
+        return str(v).rstrip('\x00')
+
     @property
     def nodata_str(self) -> str | None:
         """GDAL_NODATA tag value as string, or None."""

xrspatial/geotiff/_writer.py

@@ -51,6 +51,7 @@
     TAG_TILE_OFFSETS,
     TAG_TILE_BYTE_COUNTS,
     TAG_PREDICTOR,
+    TAG_GDAL_METADATA,
 )
 
 # Byte order: always write little-endian
@@ -412,6 +413,7 @@ def _assemble_tiff(width: int, height: int, dtype: np.dtype,
                    nodata,
                    is_cog: bool = False,
                    raster_type: int = 1,
+                   gdal_metadata_xml: str | None = None,
                    x_resolution: float | None = None,
                    y_resolution: float | None = None,
                    resolution_unit: int | None = None) -> bytes:
@@ -516,6 +518,11 @@ def _assemble_tiff(width: int, height: int, dtype: np.dtype,
                 elif gtag == TAG_GDAL_NODATA:
                     tags.append((gtag, ASCII, len(str(gval)) + 1, str(gval)))
 
+            # GDALMetadata XML (tag 42112)
+            if gdal_metadata_xml is not None:
+                tags.append((TAG_GDAL_METADATA, ASCII,
+                             len(gdal_metadata_xml) + 1, gdal_metadata_xml))
+
         ifd_specs.append(tags)
 
     # --- Determine if BigTIFF is needed ---
@@ -703,7 +710,8 @@ def write(data: np.ndarray, path: str, *,
           raster_type: int = 1,
           x_resolution: float | None = None,
           y_resolution: float | None = None,
-          resolution_unit: int | None = None) -> None:
+          resolution_unit: int | None = None,
+          gdal_metadata_xml: str | None = None) -> None:
     """Write a numpy array as a GeoTIFF or COG.
 
     Parameters
@@ -772,6 +780,7 @@ def write(data: np.ndarray, path: str, *,
         w, h, data.dtype, comp_tag, predictor, tiled, tile_size,
         parts, geo_transform, crs_epsg, nodata, is_cog=cog,
         raster_type=raster_type,
+        gdal_metadata_xml=gdal_metadata_xml,
         x_resolution=x_resolution, y_resolution=y_resolution,
         resolution_unit=resolution_unit,
     )

xrspatial/geotiff/tests/test_features.py

@@ -413,6 +413,126 @@ def test_no_crs_no_wkt(self, tmp_path):
 # Resolution / DPI tags
 # -----------------------------------------------------------------------
 
+# -----------------------------------------------------------------------
+# GDAL metadata (tag 42112)
+# -----------------------------------------------------------------------
+
+class TestGDALMetadata:
+
+    def test_parse_gdal_metadata_xml(self):
+        """XML parsing extracts dataset and per-band items."""
+        from xrspatial.geotiff._geotags import _parse_gdal_metadata
+        xml = (
+            '<GDALMetadata>\n'
+            '  <Item name="DataType">Generic</Item>\n'
+            '  <Item name="STATISTICS_MAX" sample="0">100.5</Item>\n'
+            '  <Item name="STATISTICS_MIN" sample="0">-5.2</Item>\n'
+            '  <Item name="BAND_NAME" sample="1">green</Item>\n'
+            '</GDALMetadata>\n'
+        )
+        meta = _parse_gdal_metadata(xml)
+        assert meta['DataType'] == 'Generic'
+        assert meta[('STATISTICS_MAX', 0)] == '100.5'
+        assert meta[('STATISTICS_MIN', 0)] == '-5.2'
+        assert meta[('BAND_NAME', 1)] == 'green'
+
+    def test_build_gdal_metadata_xml(self):
+        """Dict serializes back to valid XML."""
+        from xrspatial.geotiff._geotags import (
+            _build_gdal_metadata_xml, _parse_gdal_metadata)
+        meta = {
+            'DataType': 'Generic',
+            ('STATS_MAX', 0): '42.0',
+            ('STATS_MIN', 0): '-1.0',
+        }
+        xml = _build_gdal_metadata_xml(meta)
+        assert '<GDALMetadata>' in xml
+        assert '<Item name="DataType">Generic</Item>' in xml
+        assert 'sample="0"' in xml
+        # Round-trip through parser
+        reparsed = _parse_gdal_metadata(xml)
+        assert reparsed == meta
+
+    def test_round_trip_via_file(self, tmp_path):
+        """GDAL metadata survives write -> read."""
+        meta = {
+            'DataType': 'Elevation',
+            ('STATISTICS_MAXIMUM', 0): '2500.0',
+            ('STATISTICS_MINIMUM', 0): '100.0',
+            ('STATISTICS_MEAN', 0): '1200.5',
+        }
+        from xrspatial.geotiff._geotags import _build_gdal_metadata_xml
+        xml = _build_gdal_metadata_xml(meta)
+
+        arr = np.ones((4, 4), dtype=np.float32)
+        path = str(tmp_path / 'gdal_meta.tif')
+        write(arr, path, compression='none', tiled=False,
+              gdal_metadata_xml=xml)
+
+        da = read_geotiff(path)
+        assert 'gdal_metadata' in da.attrs
+        assert 'gdal_metadata_xml' in da.attrs
+        result_meta = da.attrs['gdal_metadata']
+        assert result_meta['DataType'] == 'Elevation'
+        assert result_meta[('STATISTICS_MAXIMUM', 0)] == '2500.0'
+        assert result_meta[('STATISTICS_MEAN', 0)] == '1200.5'
+
+    def test_dataarray_attrs_round_trip(self, tmp_path):
+        """GDAL metadata from DataArray attrs is preserved."""
+        meta = {'Source': 'test', ('BAND', 0): 'dem'}
+        da = xr.DataArray(
+            np.ones((4, 4), dtype=np.float32),
+            dims=['y', 'x'],
+            attrs={'gdal_metadata': meta},
+        )
+        path = str(tmp_path / 'da_meta.tif')
+        write_geotiff(da, path, compression='none')
+
+        result = read_geotiff(path)
+        assert result.attrs['gdal_metadata']['Source'] == 'test'
+        assert result.attrs['gdal_metadata'][('BAND', 0)] == 'dem'
+
+    def test_no_metadata_no_attrs(self, tmp_path):
+        """Files without GDAL metadata don't get the attrs."""
+        arr = np.ones((4, 4), dtype=np.float32)
+        path = str(tmp_path / 'no_meta.tif')
+        write(arr, path, compression='none', tiled=False)
+
+        da = read_geotiff(path)
+        assert 'gdal_metadata' not in da.attrs
+        assert 'gdal_metadata_xml' not in da.attrs
+
+    def test_real_file_metadata(self):
+        """Real USGS file has GDAL metadata with statistics."""
+        import os
+        path = '../rtxpy/examples/USGS_one_meter_x65y454_NY_LongIsland_Z18_2014.tif'
+        if not os.path.exists(path):
+            pytest.skip("Real test files not available")
+
+        da = read_geotiff(path)
+        meta = da.attrs.get('gdal_metadata')
+        assert meta is not None
+        assert 'DataType' in meta
+        assert ('STATISTICS_MAXIMUM', 0) in meta
+
+    def test_real_file_round_trip(self):
+        """GDAL metadata survives real-file round-trip."""
+        import os, tempfile
+        path = '../rtxpy/examples/USGS_one_meter_x65y454_NY_LongIsland_Z18_2014.tif'
+        if not os.path.exists(path):
+            pytest.skip("Real test files not available")
+
+        da = read_geotiff(path)
+        orig_meta = da.attrs['gdal_metadata']
+
+        out = os.path.join(tempfile.mkdtemp(), 'rt.tif')
+        write_geotiff(da, out, compression='deflate', tiled=False)
+
+        da2 = read_geotiff(out)
+        for k, v in orig_meta.items():
+            assert da2.attrs['gdal_metadata'].get(k) == v, f"Mismatch on {k}"
+
+
 class TestResolution:
 
     def test_write_read_dpi(self, tmp_path):