BUG: to_file ignores crs without any info

[raphaelquast]

• I have checked that this issue has not already been reported.

• I have confirmed this bug exists on the latest version of geopandas.

• (optional) I have confirmed this bug exists on the main branch of geopandas.

---

Code Sample, a copy-pastable example

# Your code here

from shapely.geometry import Polygon
import geopandas as gpd

# # here's an example of a crs that is silently ignored even when passed as wkt-string 
# # (e.g. the origin of the wkt string below)
# from cartopy import crs as ccrs
# crs_wkt = ccrs.PlateCarree().to_wkt()

crs_wkt = 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Equidistant Cylindrical",ID["EPSG",1028]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["unknown",111319.490793274],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["unknown",111319.490793274],ID["EPSG",8807]]],CS[Cartesian,3],AXIS["(E)",east,ORDER[1],LENGTHUNIT["unknown",111319.490793274]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["unknown",111319.490793274]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

gdf = gpd.GeoDataFrame(geometry=[Polygon([[1,2], [3,4], [5,4]])], crs=crs_wkt)

gdf.to_file("poly.shp")   # no .prj file will be created and no warning whatsoever is issued

Problem description

I guess the general problem is obvious...
The reason why I think this is problematic is because geopandas is perfectly happy with the crs (reprojections etc. work as expected) so one would expect that saving the GeoDataFrame should work as well.

Expected Output

a .proj file with the crs alongside the other shapefiles or at least a warning that the crs was not exported.

Output of geopandas.show_versions()

SYSTEM INFO

python : 3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 05:37:49) [MSC v.1916 64 bit (AMD64)]
executable : D:\miniconda\envs\rt1\python.exe
machine : Windows-10-10.0.17763-SP0

GEOS, GDAL, PROJ INFO

GEOS : None
GEOS lib : None
GDAL : 3.4.0
GDAL data dir: None
PROJ : 8.2.0
PROJ data dir: D:\miniconda\envs\rt1\Library\share\proj

PYTHON DEPENDENCIES

geopandas : 0.10.2
pandas : 1.3.5
fiona : 1.8.20
numpy : 1.21.5
shapely : 1.8.0
rtree : 0.9.7
pyproj : 3.2.1
matplotlib : 3.5.1
mapclassify: 2.4.3
geopy : None
psycopg2 : None
geoalchemy2: None
pyarrow : None
pygeos : None

[raphaelquast]

I digged a bit deeper into this and using osgeo directly, I at least get a meaningful error...
(however since its so horrifyingly undocumented, it's close to impossible to find out how to actually tell it to use "WKT2_2019" 🤯 )

from osgeo import ogr, osr

crs_wkt = 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Equidistant Cylindrical",ID["EPSG",1028]],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["unknown",111319.490793274],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["unknown",111319.490793274],ID["EPSG",8807]]],CS[Cartesian,3],AXIS["(E)",east,ORDER[1],LENGTHUNIT["unknown",111319.490793274]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["unknown",111319.490793274]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

spatialRef = osr.SpatialReference()
spatialRef.ImportFromWkt(crs_wkt)
spatialRef.MorphToESRI()
file = open(r"C:\Users\rquast\Desktop\_delete\testshp\test.prj", 'w')
file.write(spatialRef.ExportToWkt())
file.write(m.crs_plot.to_wkt())
file.close()

results in

>>> ERROR 1: PROJ: proj_as_wkt: Projected 3D CRS can only be exported since WKT2:2019

[snowman2]

This may be a helpful reference: https://pyproj4.github.io/pyproj/stable/crs_compatibility.html#osgeo-gdal

spatialRef.ExportToWkt(["FORMAT=WKT2_2018"])

[martinfleis]

Hi, thanks for the report!

I believe that this is caused by some behaviour of PROJ, that instead of raising when trying to export CRS to WKT1 GDAL format just silently returns None (see also pyproj4/pyproj#1036).

When trying to save it with pyogrio, the following warning shows up.

pyogrio.write_dataframe(gdf, 'tst.shp')

UserWarning: 'crs' was not provided.  The output dataset will not have projection information defined and may not be usable in other systems.

Even that is not entirely clear because the CRS was provided, it is just incompatible.

Because we do this conversion on our side, we should add a check and warn about that. Around here:

geopandas/geopandas/io/file.py

Lines 392 to 395 in 51864ac

	if Version(gdal_version) >= Version("3.0.0") and crs:
	crs_wkt = crs.to_wkt()
	elif crs:
	crs_wkt = crs.to_wkt("WKT1_GDAL")

[martinfleis]

This may be a helpful reference: https://pyproj4.github.io/pyproj/stable/crs_compatibility.html#osgeo-gdal

Ah, yes. Thanks!

[martinfleis]

I actually isn't that easy as with GDAL >= 3.0, we use crs.to_wkt() which has a correct output. The same situation then likely happens in fiona, which silently ignores is as it cannot be expressed as ESRI WKT (I guess).

[raphaelquast]

Hey, thanks for the follow up on this!
I just wanted to write the same thing.... my gdal-version is actually fiona.env.get_gdal_release_name() -> '3.4.0'
so WKT2 should be working just fine...

I'm really no expert on this but I think the problem might be that in pyproj there's actually only
WKT1_ESRI = 'WKT1_ESRI' pyproj docs

As a general comment, I think a warning should go to to_file in case writing the crs (e.g. a .proj file for shapefiles) fails... To me, this is the problematic part since you don't get informed that no crs-info has been stored and you just expect that things have been working fine.

[martinfleis]

As a general comment, I think a warning should go to to_file in case writing the crs (e.g. a .proj file for shapefiles) fails

Ideally yes, but right now we don't know this. We would need fiona to tell us that it has failed but it doesn't. So we would need to either guess that it is going to happen or check if .prj is present afterwards, which may be overkill.

[snowman2]

Related: Toblerity/Fiona#979

[snowman2]

Related: #1697

[snowman2]

From: Toblerity/Fiona#977 (comment)

Write:

with fiona.Env(OSR_WKT_FORMAT="WKT2_2018"), fiona.open("./test.shp", "w", driver="ESRI Shapefile", schema=schema, crs=target_crs) as collection:
    collection.write(feature):

Read:

with fiona.Env(OSR_WKT_FORMAT="WKT2_2018"), fiona.open("./test.shp") as collection:
   print(collection.crs)

[snowman2]

geopandas version: #1697 (comment)

Write:

with fiona.Env(OSR_WKT_FORMAT="WKT2_2018"):
    df.to_file('./test.shp')

Read:

with fiona.Env(OSR_WKT_FORMAT="WKT2_2018"):
   gdf = geopandas.read_file("./test.shp")