Read and write resolution/DPI tags (282, 283, 296)

Adds support for TIFF resolution metadata used in print and
cartographic workflows:

- Tag 282 (XResolution): pixels per unit, stored as RATIONAL
- Tag 283 (YResolution): pixels per unit, stored as RATIONAL
- Tag 296 (ResolutionUnit): 1=none, 2=inch, 3=centimeter

Read: resolution values are stored in DataArray attrs as
x_resolution, y_resolution (float), and resolution_unit (string:
'none', 'inch', or 'centimeter').

Write: accepted as keyword args on write() and write_geotiff(),
or extracted automatically from DataArray attrs. Written as
RATIONAL tags (numerator/denominator pairs).

Also adds RATIONAL type serialization to the writer's tag encoder.

5 new tests: DPI round-trip, centimeter unit, no-resolution check,
DataArray attrs preservation, unit='none'.

Author: brendancol

xrspatial/geotiff/__init__.py

@@ -135,6 +135,16 @@ def read_geotiff(source: str, *, window=None,
     if geo_info.raster_type == RASTER_PIXEL_IS_POINT:
         attrs['raster_type'] = 'point'
 
+    # Resolution / DPI metadata
+    if geo_info.x_resolution is not None:
+        attrs['x_resolution'] = geo_info.x_resolution
+    if geo_info.y_resolution is not None:
+        attrs['y_resolution'] = geo_info.y_resolution
+    if geo_info.resolution_unit is not None:
+        _unit_names = {1: 'none', 2: 'inch', 3: 'centimeter'}
+        attrs['resolution_unit'] = _unit_names.get(
+            geo_info.resolution_unit, str(geo_info.resolution_unit))
+
     # Attach palette colormap for indexed-color TIFFs
     if geo_info.colormap is not None:
         try:
@@ -219,6 +229,9 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
     geo_transform = None
     epsg = crs
     raster_type = RASTER_PIXEL_IS_AREA
+    x_res = None
+    y_res = None
+    res_unit = None
 
     if isinstance(data, xr.DataArray):
         arr = data.values
@@ -230,6 +243,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
+        # Resolution / DPI from attrs
+        x_res = data.attrs.get('x_resolution')
+        y_res = data.attrs.get('y_resolution')
+        unit_str = data.attrs.get('resolution_unit')
+        if unit_str is not None:
+            _unit_ids = {'none': 1, 'inch': 2, 'centimeter': 3}
+            res_unit = _unit_ids.get(str(unit_str), None)
     else:
         arr = np.asarray(data)
 
@@ -249,6 +269,9 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
         overview_levels=overview_levels,
         overview_resampling=overview_resampling,
         raster_type=raster_type,
+        x_resolution=x_res,
+        y_resolution=y_res,
+        resolution_unit=res_unit,
     )
 
 

xrspatial/geotiff/_geotags.py

@@ -63,6 +63,9 @@ class GeoInfo:
     raster_type: int = RASTER_PIXEL_IS_AREA
     nodata: float | None = None
     colormap: list | None = None  # list of (R, G, B, A) float tuples, or None
+    x_resolution: float | None = None
+    y_resolution: float | None = None
+    resolution_unit: int | None = None  # 1=none, 2=inch, 3=cm
     geokeys: dict[int, int | float | str] = field(default_factory=dict)
 
 
@@ -267,6 +270,9 @@ def extract_geo_info(ifd: IFD, data: bytes | memoryview,
         raster_type=int(raster_type) if isinstance(raster_type, (int, float)) else RASTER_PIXEL_IS_AREA,
         nodata=nodata,
         colormap=colormap,
+        x_resolution=ifd.x_resolution,
+        y_resolution=ifd.y_resolution,
+        resolution_unit=ifd.resolution_unit,
         geokeys=geokeys,
     )
 

xrspatial/geotiff/_header.py

@@ -24,7 +24,10 @@
 TAG_SAMPLES_PER_PIXEL = 277
 TAG_ROWS_PER_STRIP = 278
 TAG_STRIP_BYTE_COUNTS = 279
+TAG_X_RESOLUTION = 282
+TAG_Y_RESOLUTION = 283
 TAG_PLANAR_CONFIG = 284
+TAG_RESOLUTION_UNIT = 296
 TAG_PREDICTOR = 317
 TAG_TILE_WIDTH = 322
 TAG_TILE_LENGTH = 323
@@ -159,6 +162,23 @@ def photometric(self) -> int:
     def planar_config(self) -> int:
         return self.get_value(TAG_PLANAR_CONFIG, 1)
 
+    @property
+    def x_resolution(self) -> float | None:
+        """XResolution tag (282), or None if absent."""
+        v = self.get_value(TAG_X_RESOLUTION)
+        return float(v) if v is not None else None
+
+    @property
+    def y_resolution(self) -> float | None:
+        """YResolution tag (283), or None if absent."""
+        v = self.get_value(TAG_Y_RESOLUTION)
+        return float(v) if v is not None else None
+
+    @property
+    def resolution_unit(self) -> int | None:
+        """ResolutionUnit tag (296): 1=none, 2=inch, 3=cm. None if absent."""
+        return self.get_value(TAG_RESOLUTION_UNIT)
+
     @property
     def colormap(self) -> tuple | None:
         """ColorMap tag (320) values, or None if absent."""

xrspatial/geotiff/_writer.py

@@ -17,6 +17,7 @@
 )
 from ._dtypes import (
     DOUBLE,
+    RATIONAL,
     SHORT,
     LONG,
     ASCII,
@@ -42,6 +43,9 @@
     TAG_STRIP_OFFSETS,
     TAG_ROWS_PER_STRIP,
     TAG_STRIP_BYTE_COUNTS,
+    TAG_X_RESOLUTION,
+    TAG_Y_RESOLUTION,
+    TAG_RESOLUTION_UNIT,
     TAG_TILE_WIDTH,
     TAG_TILE_LENGTH,
     TAG_TILE_OFFSETS,
@@ -154,6 +158,16 @@ def _make_overview(arr: np.ndarray, method: str = 'mean') -> np.ndarray:
 # Tag serialization
 # ---------------------------------------------------------------------------
 
+def _float_to_rational(val):
+    """Convert a float to a TIFF RATIONAL (numerator, denominator) pair."""
+    if val == int(val):
+        return (int(val), 1)
+    # Use a denominator of 10000 for reasonable precision
+    den = 10000
+    num = int(round(val * den))
+    return (num, den)
+
+
 def _serialize_tag_value(type_id, count, values):
     """Serialize tag values to bytes."""
     if type_id == ASCII:
@@ -168,6 +182,16 @@ def _serialize_tag_value(type_id, count, values):
         if isinstance(values, (list, tuple)):
             return struct.pack(f'{BO}{count}I', *values)
         return struct.pack(f'{BO}I', values)
+    elif type_id == RATIONAL:
+        # RATIONAL = two LONGs (numerator, denominator) per value
+        if isinstance(values, (list, tuple)) and isinstance(values[0], (list, tuple)):
+            parts = []
+            for num, den in values:
+                parts.extend([int(num), int(den)])
+            return struct.pack(f'{BO}{count * 2}I', *parts)
+        else:
+            num, den = _float_to_rational(float(values))
+            return struct.pack(f'{BO}II', num, den)
     elif type_id == DOUBLE:
         if isinstance(values, (list, tuple)):
             return struct.pack(f'{BO}{count}d', *values)
@@ -387,7 +411,10 @@ def _assemble_tiff(width: int, height: int, dtype: np.dtype,
                    crs_epsg: int | None,
                    nodata,
                    is_cog: bool = False,
-                   raster_type: int = 1) -> bytes:
+                   raster_type: int = 1,
+                   x_resolution: float | None = None,
+                   y_resolution: float | None = None,
+                   resolution_unit: int | None = None) -> bytes:
     """Assemble a complete TIFF file.
 
     Parameters
@@ -455,6 +482,14 @@ def _assemble_tiff(width: int, height: int, dtype: np.dtype,
         if pred_val != 1:
             tags.append((TAG_PREDICTOR, SHORT, 1, pred_val))
 
+        # Resolution / DPI tags
+        if x_resolution is not None:
+            tags.append((TAG_X_RESOLUTION, RATIONAL, 1, x_resolution))
+        if y_resolution is not None:
+            tags.append((TAG_Y_RESOLUTION, RATIONAL, 1, y_resolution))
+        if resolution_unit is not None:
+            tags.append((TAG_RESOLUTION_UNIT, SHORT, 1, resolution_unit))
+
         if tiled:
             tags.append((TAG_TILE_WIDTH, SHORT, 1, tile_size))
             tags.append((TAG_TILE_LENGTH, SHORT, 1, tile_size))
@@ -665,7 +700,10 @@ def write(data: np.ndarray, path: str, *,
           cog: bool = False,
           overview_levels: list[int] | None = None,
           overview_resampling: str = 'mean',
-          raster_type: int = 1) -> None:
+          raster_type: int = 1,
+          x_resolution: float | None = None,
+          y_resolution: float | None = None,
+          resolution_unit: int | None = None) -> None:
     """Write a numpy array as a GeoTIFF or COG.
 
     Parameters
@@ -734,6 +772,8 @@ 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,
+        x_resolution=x_resolution, y_resolution=y_resolution,
+        resolution_unit=resolution_unit,
     )
 
     # Write to a temp file then atomically rename, so concurrent writes to

xrspatial/geotiff/tests/test_features.py

@@ -256,6 +256,73 @@ def test_zstd_public_api(self, tmp_path):
         np.testing.assert_array_equal(result.values, arr)
 
 
+# -----------------------------------------------------------------------
+# Resolution / DPI tags
+# -----------------------------------------------------------------------
+
+class TestResolution:
+
+    def test_write_read_dpi(self, tmp_path):
+        """Resolution tags round-trip through write and read."""
+        arr = np.ones((4, 4), dtype=np.float32)
+        path = str(tmp_path / 'dpi.tif')
+        write(arr, path, compression='none', tiled=False,
+              x_resolution=300.0, y_resolution=300.0, resolution_unit=2)
+
+        da = read_geotiff(path)
+        assert da.attrs['x_resolution'] == pytest.approx(300.0, rel=0.01)
+        assert da.attrs['y_resolution'] == pytest.approx(300.0, rel=0.01)
+        assert da.attrs['resolution_unit'] == 'inch'
+
+    def test_write_read_cm(self, tmp_path):
+        """Centimeter resolution unit."""
+        arr = np.ones((4, 4), dtype=np.float32)
+        path = str(tmp_path / 'dpi_cm.tif')
+        write(arr, path, compression='none', tiled=False,
+              x_resolution=118.0, y_resolution=118.0, resolution_unit=3)
+
+        da = read_geotiff(path)
+        assert da.attrs['x_resolution'] == pytest.approx(118.0, rel=0.01)
+        assert da.attrs['resolution_unit'] == 'centimeter'
+
+    def test_no_resolution_no_attrs(self, tmp_path):
+        """Files without resolution tags don't get resolution attrs."""
+        arr = np.ones((4, 4), dtype=np.float32)
+        path = str(tmp_path / 'no_dpi.tif')
+        write(arr, path, compression='none', tiled=False)
+
+        da = read_geotiff(path)
+        assert 'x_resolution' not in da.attrs
+        assert 'y_resolution' not in da.attrs
+        assert 'resolution_unit' not in da.attrs
+
+    def test_dataarray_attrs_round_trip(self, tmp_path):
+        """Resolution attrs on DataArray are preserved through write/read."""
+        da = xr.DataArray(
+            np.ones((4, 4), dtype=np.float32),
+            dims=['y', 'x'],
+            attrs={'x_resolution': 72.0, 'y_resolution': 72.0,
+                   'resolution_unit': 'inch'},
+        )
+        path = str(tmp_path / 'da_dpi.tif')
+        write_geotiff(da, path, compression='none')
+
+        result = read_geotiff(path)
+        assert result.attrs['x_resolution'] == pytest.approx(72.0, rel=0.01)
+        assert result.attrs['y_resolution'] == pytest.approx(72.0, rel=0.01)
+        assert result.attrs['resolution_unit'] == 'inch'
+
+    def test_unit_none(self, tmp_path):
+        """ResolutionUnit=1 (no unit) round-trips as 'none'."""
+        arr = np.ones((4, 4), dtype=np.float32)
+        path = str(tmp_path / 'no_unit.tif')
+        write(arr, path, compression='none', tiled=False,
+              x_resolution=1.0, y_resolution=1.0, resolution_unit=1)
+
+        da = read_geotiff(path)
+        assert da.attrs['resolution_unit'] == 'none'
+
+
 # -----------------------------------------------------------------------
 # Overview resampling methods
 # -----------------------------------------------------------------------