Add support for reading and writing tiff gcps

rioxarray/_io.py

@@ -533,9 +533,37 @@ def _get_rasterio_attrs(riods):
         else:
             attrs["units"] = riods.units
 
+    if hasattr(riods, "gcps") and riods.gcps != ([], None):
+        gcps, gcp_crs = riods.gcps
+        attrs["gcps"] = _convert_gcps_to_geojson(gcps)
+        attrs["gcp_crs"] = gcp_crs
+
     return attrs
 
 
+def _convert_gcps_to_geojson(gcps):
+    """
+    Convert GCPs to geojson.
+
+    Parameters
+    ----------
+    gcps: The list of GroundControlPoint instances.
+
+    Returns
+    -------
+    A FeatureCollection dict.
+    """
+    features = [
+        dict(
+            type="Feature",
+            properties=dict(id=gcp.id, info=gcp.info, row=gcp.row, col=gcp.col),
+            geometry=dict(type="Point", coordinates=[gcp.x, gcp.y, gcp.z]),
+        )
+        for gcp in gcps
+    ]
+    return dict(type="FeatureCollection", features=features)
+
+
 def _decode_datetime_cf(data_array, decode_times, decode_timedelta):
     """
     Decide the datetime based on CF conventions

test/integration/test_integration__io.py

@@ -760,7 +760,7 @@ def test_ENVI_tags():
 
 
 def test_no_mftime():
-    # rasterio can accept "filename" urguments that are actually urls,
+    # rasterio can accept "filename" arguments that are actually urls,
     # including paths to remote files.
     # In issue #1816, we found that these caused dask to break, because
     # the modification time was used to determine the dask token. This
@@ -1236,3 +1236,52 @@ def test_cint16_promote_dtype(tmp_path):
         assert "complex128" in riofh.dtypes
         assert riofh.nodata == 0
         assert data.dtype == "complex128"
+
+
+def test_reading_gcps(tmp_path):
+    """
+    Test reading gcps from a tiff file.
+    """
+    tiffname = tmp_path / "test.tif"
+    src_gcps = [
+        GroundControlPoint(
+            row=0, col=0, x=0.0, y=0.0, z=12.0, id="1", info="the first gcp"
+        ),
+        GroundControlPoint(
+            row=0, col=800, x=10.0, y=0.0, z=1.0, id="2", info="the second gcp"
+        ),
+        GroundControlPoint(
+            row=800, col=800, x=10.0, y=10.0, z=3.5, id="3", info="the third gcp"
+        ),
+        GroundControlPoint(
+            row=800, col=0, x=0.0, y=10.0, z=5.5, id="4", info="the fourth gcp"
+        ),
+    ]
+    crs = CRS.from_epsg(4326)
+    with rasterio.open(
+        tiffname,
+        mode="w",
+        height=800,
+        width=800,
+        count=3,
+        dtype=np.uint8,
+        driver="GTiff",
+    ) as source:
+        source.gcps = (src_gcps, crs)
+
+    with rioxarray.open_rasterio(tiffname) as darr:
+        assert "gcps" in darr.attrs
+        assert "gcp_crs" in darr.attrs
+        assert darr.attrs["gcp_crs"] == crs
+        gcps = darr.attrs["gcps"]
+        assert gcps["type"] == "FeatureCollection"
+        assert len(gcps["features"]) == len(src_gcps)
+        for feature, gcp in zip(gcps["features"], src_gcps):
+            assert feature["type"] == "Feature"
+            assert feature["properties"]["id"] == gcp.id
+            # info seems to be lost when rasterio writes?
+            # assert feature["properties"]["info"] == gcp.info
+            assert feature["properties"]["row"] == gcp.row
+            assert feature["properties"]["col"] == gcp.col
+            assert feature["geometry"]["type"] == "Point"
+            assert feature["geometry"]["coordinates"] == [gcp.x, gcp.y, gcp.z]