2

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]]
  • 1
    Which python library or tools do you use to read your shapefile? – YoLecomte Nov 27 '18 at 11:21
  • I have tried GeoPandas, Shapely, and Fionna – userxxx Nov 27 '18 at 11:34
3

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

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

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

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

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

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.