How to Identify Coordinate System (CRS) from projection file using Python

Question

I am working on many shapefiles (their corresponding projection files are also presented).

Every shapefile has a different Coordinate System, so it becomes difficult to plot these in a generic manner.

How can I identify with Python what is the CRS used in the shapefile?

For example, the ".prj" file for some shapefiles is given below, but I am not able to identify what is the CRS to be used.

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

GEOGCS["WGS84(DD)", DATUM["WGS84", SPHEROID["WGS84", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH]]

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

GEOGCS["WGS84(DD)", DATUM["WGS84", SPHEROID["WGS84", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH]]

Which python library or tools do you use to read your shapefile? — YoLecomte, Nov 27 '18 at 11:21

YoLecomte · Answer 1 · 2018-11-27T11:52:50.727

6

For example with fiona:

c = fiona.open('docs/data/test_uk.shp')
crs = c.crs

find it in the docs page 28

If you need the EPSG code there is good inspiration to take here using pyproj

edited Nov 27 '18 at 11:52

answered Nov 27 '18 at 11:44

YoLecomte

2,895
10
23

score 2 · Answer 2 · answered Nov 27 '18 at 17:10

A number of Python modules exist to work with CRS or projection files

Pypi:ESPSG or Pypi:CRS
GitHub: EPSG

For example, you can use epsg-ident or sridentify (Quickly get the EPSG code from a .prj file or WKT)

from sridentify import Sridentify
ident = Sridentify()
# from file
ident.from_file('schisto.prj')
ident.get_epsg()
31370
# from WKT
ident = Sridentify(prj="""GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]])
ident.get_epsg()
4326

But it doesn't always work

ident =  Sridentify(prj = """GEOGCS["WGS84(DD)", DATUM["WGS84", SPHEROID["WGS84", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH]""")
ident.get_epsg()
# nothing

score 1 · Answer 3 · answered Nov 27 '18 at 11:28

1

Is it not just that simple as below (I know he wanted to know from prj file but I think this is more simple)?

import arcpy

arcpy.env.workspace = r"PATH TO FOLDER OR GDB"

featureclasses = arcpy.ListFeatureClasses()
for fc in featureclasses:
    desc = arcpy.Describe(fc)
    spatialRef = desc.spatialReference
    print(spatialRef.Name)

https://www.e-education.psu.edu/geog485/node/115

answered Nov 27 '18 at 11:28

JuniorPythonNewbie

573
2
12

Thanks for this, but is there any alternate solution which does not use arcpy. – userxxx Nov 27 '18 at 11:37

score 0 · Answer 4 · answered Jan 20 '22 at 12:01

0

In two lines using geopandas (pip install geopandas):

import geopandas as gpd
print(gpd.read_file('world-administrative-boundaries.shp').crs)

answered Jan 20 '22 at 12:01

G M

2,068
2
21
31

How to Identify Coordinate System (CRS) from projection file using Python

4 Answers4