ENH: include crs in to_json (#1774) (#2151)

Co-authored-by: Jan Šimbera <jan.simbera@nanoenergies.cz>
Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>
Co-authored-by: Matt Richards <45483497+m-richards@users.noreply.github.com>
Co-authored-by: Joris Van den Bossche <jorisvandenbossche@gmail.com>

CHANGELOG.md

@@ -19,6 +19,7 @@ New features and improvements:
   using ``engine="pyogrio"`` (#2788).
 - Added a ``to_wgs84`` keyword to ``to_json`` allowing automatic re-projecting to follow
   the 2016 GeoJSON specification (#416).
+- ``to_json`` output now includes a ``"crs"`` field if the CRS is not the default WGS84 (#1774).
 - Improve error messages when accessing the `geometry` attribute of GeoDataFrame without an active geometry column 
   related to the default name `"geometry"` being provided in the constructor (#2577)
 

geopandas/geodataframe.py

@@ -766,22 +766,28 @@ def to_json(
           feature individually so that features may have different properties.
         - ``keep``: output the missing entries as NaN.
 
+        If the GeoDataFrame has a defined CRS, its definition will be included
+        in the output unless it is equal to WGS84 (default GeoJSON CRS) or not
+        possible to represent in the URN OGC format, or unless ``to_wgs84=True``
+        is specified.
+
         Examples
         --------
 
         >>> from shapely.geometry import Point
         >>> d = {'col1': ['name1', 'name2'], 'geometry': [Point(1, 2), Point(2, 1)]}
-        >>> gdf = geopandas.GeoDataFrame(d, crs="EPSG:4326")
+        >>> gdf = geopandas.GeoDataFrame(d, crs="EPSG:3857")
         >>> gdf
-            col1                 geometry
-        0  name1  POINT (1.00000 2.00000)
-        1  name2  POINT (2.00000 1.00000)
+            col1             geometry
+        0  name1  POINT (1.000 2.000)
+        1  name2  POINT (2.000 1.000)
 
         >>> gdf.to_json()
         '{"type": "FeatureCollection", "features": [{"id": "0", "type": "Feature", \
 "properties": {"col1": "name1"}, "geometry": {"type": "Point", "coordinates": [1.0,\
  2.0]}}, {"id": "1", "type": "Feature", "properties": {"col1": "name2"}, "geometry"\
-: {"type": "Point", "coordinates": [2.0, 1.0]}}]}'
+: {"type": "Point", "coordinates": [2.0, 1.0]}}], "crs": {"type": "name", "properti\
+es": {"name": "urn:ogc:def:crs:EPSG::3857"}}}'
 
         Alternatively, you can write GeoJSON to file:
 
@@ -801,9 +807,26 @@ def to_json(
                 )
         else:
             df = self
-        return json.dumps(
-            df._to_geo(na=na, show_bbox=show_bbox, drop_id=drop_id), **kwargs
-        )
+
+        geo = df._to_geo(na=na, show_bbox=show_bbox, drop_id=drop_id)
+
+        # if the geometry is not in WGS84, include CRS in the JSON
+        if df.crs is not None and not df.crs.equals("epsg:4326"):
+            auth_crsdef = self.crs.to_authority()
+            allowed_authorities = ["EDCS", "EPSG", "OGC", "SI", "UCUM"]
+
+            if auth_crsdef is None or auth_crsdef[0] not in allowed_authorities:
+                warnings.warn(
+                    "GeoDataFrame's CRS is not representable in URN OGC "
+                    "format. Resulting JSON will contain no CRS information.",
+                    stacklevel=2,
+                )
+            else:
+                authority, code = auth_crsdef
+                ogc_crs = f"urn:ogc:def:crs:{authority}::{code}"
+                geo["crs"] = {"type": "name", "properties": {"name": ogc_crs}}
+
+        return json.dumps(geo, **kwargs)
 
     @property
     def __geo_interface__(self):
@@ -814,7 +837,10 @@ def __geo_interface__(self):
         ``FeatureCollection``.
 
         This differs from `_to_geo()` only in that it is a property with
-        default args instead of a method
+        default args instead of a method.
+
+        CRS of the dataframe is not passed on to the output, unlike
+        :meth:`~GeoDataFrame.to_json()`.
 
         Examples
         --------
@@ -833,8 +859,6 @@ def __geo_interface__(self):
 , 2.0)}, 'bbox': (1.0, 2.0, 1.0, 2.0)}, {'id': '1', 'type': 'Feature', 'properties\
 ': {'col1': 'name2'}, 'geometry': {'type': 'Point', 'coordinates': (2.0, 1.0)}, 'b\
 box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
-
-
         """
         return self._to_geo(na="null", show_bbox=True, drop_id=False)
 

geopandas/tests/test_geodataframe.py

@@ -1427,3 +1427,45 @@ def test_geodataframe_crs():
     gdf = GeoDataFrame(columns=["geometry"])
     gdf.crs = "IGNF:ETRS89UTM28"
     assert gdf.crs.to_authority() == ("IGNF", "ETRS89UTM28")
+
+
+def test_geodataframe_nocrs_json():
+    # no CRS, no crs field
+    gdf = GeoDataFrame(columns=["geometry"])
+    gdf_geojson = json.loads(gdf.to_json())
+    assert "crs" not in gdf_geojson
+
+    # WGS84, no crs field (default as per spec)
+    gdf.crs = 4326
+    gdf_geojson = json.loads(gdf.to_json())
+    assert "crs" not in gdf_geojson
+
+
+def test_geodataframe_crs_json():
+    gdf = GeoDataFrame(columns=["geometry"])
+    gdf.crs = 25833
+    gdf_geojson = json.loads(gdf.to_json())
+    assert "crs" in gdf_geojson
+    assert gdf_geojson["crs"] == {
+        "type": "name",
+        "properties": {"name": "urn:ogc:def:crs:EPSG::25833"},
+    }
+    gdf_geointerface = gdf.__geo_interface__
+    assert "crs" not in gdf_geointerface
+
+
+@pytest.mark.parametrize(
+    "crs",
+    ["+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m", "IGNF:WGS84"],
+)
+def test_geodataframe_crs_nonrepresentable_json(crs):
+    gdf = GeoDataFrame(
+        [Point(1000, 1000)],
+        columns=["geometry"],
+        crs=crs,
+    )
+    with pytest.warns(
+        UserWarning, match="GeoDataFrame's CRS is not representable in URN OGC"
+    ):
+        gdf_geojson = json.loads(gdf.to_json())
+    assert "crs" not in gdf_geojson