Repository: corteva/rioxarray

Issue:
Having dict-encoded GCPs leads to `to_netcdf` failing with a `TypeError`
#### Code Sample

```python
import xarray as xr
import rioxarray
from rasterio.control import GroundControlPoint
from rasterio.crs import CRS

ds = xr.Dataset(dict(A=xr.DataArray([1.0])))

gcps = [GroundControlPoint(row=0, col=0, x=10, y=10, z=0, id='0', info='')]
crs = CRS.from_epsg(4326)

ds.rio.write_gcps(gcps, crs, inplace=True)

ds.to_netcdf("/path/to/dataset.nc")
# TypeError: Invalid value for attr 'gcps': {'type': 'FeatureCollection', 'features': [{'type': 'Feature', 'properties': {'id': '0', 'info': '', 'row': 0, 'col': 0}, 'geometry': {'type': 'Point', 'coordinates': [10, 10, 0]}}]}. For serialization to netCDF files, its value must be of one of the following types: str, Number, ndarray, number, list, tuple
```

#### Problem description

Since rioxarray currently encodes GCPs as a GeoJSON dict in the `spatial_ref` attribute `gcps`, a GCP referenced dataset can not be saved as NetCDF.

#### Expected Output

I would expect to be able to save a GCP referenced dataset as NetCDF. This could easily be fixed by calling `json.dumps` on the GeoJSON-dict and storing the resulting string in the `gcps` attribute instead (and calling `json.loads` before reading the GCPs in the accessor function `get_gcps`).

Discussion:

Patch:
diff --git a/docs/history.rst b/docs/history.rst
index 7819e42e..40237617 100644
--- a/docs/history.rst
+++ b/docs/history.rst
@@ -4,6 +4,7 @@ History
 Latest
 ------
 - BUG: Raise OverflowError when nodata data type conversion is unsafe (pull #782)
+- BUG: Support writing GCPs to netCDF (issue #778)
 
 0.15.5
 ------
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 4c69783f..66fd3eac 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -3,6 +3,7 @@ This module is an extension for xarray to provide rasterio capabilities
 to xarray datasets/dataarrays.
 """
 # pylint: disable=too-many-lines
+import json
 import math
 import warnings
 from collections.abc import Hashable, Iterable
@@ -254,6 +255,7 @@ def _order_bounds(
 class XRasterBase:
     """This is the base class for the GIS extensions for xarray"""
 
+    # pylint: disable=too-many-instance-attributes
     def __init__(self, xarray_obj: Union[xarray.DataArray, xarray.Dataset]):
         self._obj: Union[xarray.DataArray, xarray.Dataset] = xarray_obj
 
@@ -289,6 +291,7 @@ class XRasterBase:
         self._height: Optional[int] = None
         self._width: Optional[int] = None
         self._crs: Union[rasterio.crs.CRS, None, Literal[False]] = None
+        self._gcps: Optional[list[GroundControlPoint]] = None
 
     @property
     def crs(self) -> Optional[rasterio.crs.CRS]:
@@ -347,6 +350,7 @@ class XRasterBase:
         obj_copy.rio._width = self._width
         obj_copy.rio._height = self._height
         obj_copy.rio._crs = self._crs
+        obj_copy.rio._gcps = self._gcps
         return obj_copy
 
     def set_crs(
@@ -1235,7 +1239,8 @@ class XRasterBase:
             gcp_crs, grid_mapping_name=grid_mapping_name, inplace=inplace
         )
         geojson_gcps = _convert_gcps_to_geojson(gcps)
-        data_obj.coords[grid_mapping_name].attrs["gcps"] = geojson_gcps
+        data_obj.coords[grid_mapping_name].attrs["gcps"] = json.dumps(geojson_gcps)
+        self._gcps = list(gcps)
         return data_obj
 
     def get_gcps(self) -> Optional[list[GroundControlPoint]]:
@@ -1249,8 +1254,10 @@ class XRasterBase:
         list of :obj:`rasterio.control.GroundControlPoint` or None
             The Ground Control Points from the dataset or None if not applicable
         """
+        if self._gcps is not None:
+            return self._gcps
         try:
-            geojson_gcps = self._obj.coords[self.grid_mapping].attrs["gcps"]
+            geojson_gcps = json.loads(self._obj.coords[self.grid_mapping].attrs["gcps"])
         except (KeyError, AttributeError):
             return None
 
@@ -1267,8 +1274,8 @@ class XRasterBase:
                 info=gcp["properties"]["info"],
             )
 
-        gcps = [_parse_gcp(gcp) for gcp in geojson_gcps["features"]]
-        return gcps
+        self._gcps = [_parse_gcp(gcp) for gcp in geojson_gcps["features"]]
+        return self._gcps
 
 
 def _convert_gcps_to_geojson(


Test patch:
diff --git a/test/integration/test_integration__io.py b/test/integration/test_integration__io.py
index b37df614..9fd4e3e4 100644
--- a/test/integration/test_integration__io.py
+++ b/test/integration/test_integration__io.py
@@ -1459,6 +1459,34 @@ def test_writing_gcps(tmp_path):
         _check_rio_gcps(darr, *gdal_gcps)
 
 
+def test_writing_gcps__to_netcdf(tmp_path):
+    """
+    Test writing gcps to a netCDF file.
+    """
+    tiffname = tmp_path / "test.tif"
+    nc_name = tmp_path / "test_written.nc"
+
+    src_gcps, crs = _create_gdal_gcps()
+
+    with rasterio.open(
+        tiffname,
+        mode="w",
+        height=800,
+        width=800,
+        count=3,
+        dtype=numpy.uint8,
+        driver="GTiff",
+    ) as source:
+        source.gcps = (src_gcps, crs)
+
+    with rioxarray.open_rasterio(tiffname) as darr:
+        darr.to_netcdf(nc_name)
+
+    with xarray.open_dataset(nc_name, decode_coords="all") as darr:
+        assert "gcps" in darr.coords["spatial_ref"].attrs
+        _check_rio_gcps(darr, src_gcps=src_gcps, crs=crs)
+
+
 def test_read_file_handle_with_dask():
     with open(
         os.path.join(TEST_COMPARE_DATA_DIR, "small_dem_3m_merged.tif"), "rb"
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index fdb3f6bf..ded645b7 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -2987,7 +2987,7 @@ def _check_rio_gcps(darr, src_gcps, crs):
     assert "y" not in darr.coords
     assert darr.rio.crs == crs
     assert "gcps" in darr.spatial_ref.attrs
-    gcps = darr.spatial_ref.attrs["gcps"]
+    gcps = json.loads(darr.spatial_ref.attrs["gcps"])
     assert gcps["type"] == "FeatureCollection"
     assert len(gcps["features"]) == len(src_gcps)
     for feature, gcp in zip(gcps["features"], src_gcps):