Repository: geopandas/geopandas

Issue:
ENH: include CRS in "to_json"
_Note (by @martinfleis): see the discussion below for details on the enhancement proposal._

#### Question about geopandas

gdf.to_json() 

> {"type": "FeatureCollection", 
> "features": [{**"id": "0500"**, "type": "Feature", **"properties": {}**, "geometry": ...}]
> }
> 

gdf.to_file('temp.geojson', driver='GeoJSON')

> {
> "type": "FeatureCollection",
> **"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32246" } },**
> "features": [
> { "type": "Feature", **"properties": { "id": "0500" }**, "geometry": ...}]
> }


the   "id" is in "properties" for .to_file(driver='GeoJSON')  
and to_json miss the "crs"

Why they have different behavior?

and

How to:

1  let to_json()  output "crs" ?

2 let to_file()  output   each feature  with "id" as  to_json()?

Thanks.

Discussion:
martinfleis: Thanks!

I am only worried about one thing. `__geo_interface__` spec does not include crs. If we now include it, following the GeoJSON model, I am not sure that it doesn't cause some trouble in other libraries. They may not expect it and hence not know how to handle it.

I would be more comfortable if we included it only in `to_json` where it is expected and not in `__geo_interface__` where it is not. 
simberaj: @martinfleis Updated according to your idea.
simberaj: @martinfleis docstring amended, thanks for the suggestion
martinfleis: I just realised that this fails when a CRS cannot be represented by an (authority, code) tuple. Consider this example, using Gall-Peters Orthographic Projection defined via proj string. The CRS is perfectly valid, just pyproj cannot determine authority and code.

```py
import geopandas

world = geopandas.read_file(geopandas.datasets.get_path('nybb'))
world = world.set_geometry(world.centroid)
world = world.set_crs('+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs', allow_override=True)
world.to_json()
```

```py
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/var/folders/2f/fhks6w_d0k556plcv3rfmshw0000gn/T/ipykernel_24983/3825038059.py in <module>
----> 1 world.to_json()

~/Git/geopandas/geopandas/geodataframe.py in to_json(self, na, show_bbox, drop_id, **kwargs)
    749 
    750         if self.crs is not None:
--> 751             authority, code = self.crs.to_authority()
    752             ogc_crs = f"urn:ogc:def:crs:{authority}::{code}"
    753             geo["crs"] = {"type": "name", "properties": {"name": ogc_crs}}

TypeError: cannot unpack non-iterable NoneType object
```

We should catch this, raise a warning that we were unable to create a compliant CRS string and return a JSON without CRS. 
martinfleis: Thanks for the change! I have pushed two commits in, as I initially just wanted to lint the code but then realised that we should remove the print of the CRS from the error message as it gets quite verbose, so I went ahead. We don't print the whole thing anywhere and I think we can easily leave it out here.
jorisvandenbossche: I don't know how important this is, but one potential additional concern: according to https://www.ogc.org/ogcUrnPolicy, there is a specific set of authorities that are currently recognized by the OGC spec (EDCS, EPSG, OGC, SI, UCUM, so eg not IGNF or ESRI). 
(might be worth checking what GDAL does in those cases)
jorisvandenbossche: Quickly testing:

```
df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
df2 = df.to_crs("IGNF:WGS84")
df2.to_file("test_ignf.geojson")
```

this then doesn't seem to include the CRS information in the resulting file (while it does for `df.to_file("test.geojson")`).
martinfleis: I have updated this PR to closely match the 2016 GeoJSON spec. That means that when we save a gdf with a CRS equal to WGS84 or use `to_wgs84=True`, no CRS is saved and the JSON matches the spec.

When we use different CRS, a `crs` field is included matching the older spec and GDAL behaviour. 

We also check for the set of allowed authorities and include CRS only if the authority is listed in the [OGC URN Policy](https://www.ogc.org/about-ogc/policies/ogc-urn-policy/). This **deviates from GDAL** which saves CRS using the authority::code even for authorities like IGNF or ESRI. But it follows OGC. The question here is whether we want to do this - the resulting GeoJSON is not valid per the specification anyway (different CRS used and no CRS stored), so we may want to follow GDAL in saving the CRS anytime we are able to obtain authority:code tuple from pyproj. 

Removing the check for `allowed_authorities`, the behaviour should match the one of GDAL, based on my testing.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 93aef7b3..1dfa6a61 100644
--- a/CHANGELOG.md
+++ b/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)
 
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 3557a2fd..3fb64d08 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -766,22 +766,28 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
           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 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
                 )
         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 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
         ``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 @@ class GeoDataFrame(GeoPandasBase, DataFrame):
 , 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)
 


Test patch:
diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py
index 9e13af72..34412706 100644
--- a/geopandas/tests/test_geodataframe.py
+++ b/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