1

Loving fiona, but have a small coordinate system issue.

I'm trying to read the crs property on a shapefile to validate projection. Below I contrast fiona documentation with what is output on my console.

Either output would be fine for my use case, but any idea why they are different?

Documenation states: http://toblerity.org/fiona/manual.html#format-drivers-crs-bounds-and-schema

>>> c = fiona.open('docs/data/test_uk.shp')
>>> c.crs
>>> {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'longlat'}

I get, as if I were using from_epsg:

>>> {'init': u'epsg:4326'}

Are there scenarios where different shapefiles / prj formats yield different crs objects?

My .prj file: GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Edit Addition: In looking at an additional shapefile with prj: PROJCS["USA_Contiguous_Albers_Equal_Area_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-96.0],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2",45.5],PARAMETER["Latitude_Of_Origin",37.5],UNIT["Meter",1.0]]

I get the following output from c.crs:

{u'lon_0': -96, u'datum': u'NAD83', u'y_0': 0, u'no_defs': True, u'proj': u'aea', u'x_0': 0, u'units': u'm', u'lat_2': 45.5, u'lat_1': 29.5, u'lat_0': 37.5}

If I am just validating projection, is there a uniform interface besides .crs to get that info? Thanks for any help

  • The crs format returned by Fiona is the Proj.4 format not the ESRI prj format – gene Nov 1 '14 at 21:42
1

Just a difference between documentation and code. Recent versions of Fiona sniff out whether there's an entry in the EPSG table corresponding to the projection and preferentially return that.

  • 1
    And in Fiona 1.4.8, feature collections have a crs_wkt property that gives WKT strings. – sgillies Nov 2 '14 at 14:01
0

Not sure, it's better to use Fiona here.

You can compare strings from your prj file for each layer using http://prj2epsg.org/search

You can also try to use OGR from Python like below (except last line, it's an extraction from the Python GDAL/OGR Cookbook

from osgeo import ogr, osr
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open(r'c:\data\yourshpfile.shp')

# from Layer
layer = dataset.GetLayer()
spatialRef = layer.GetSpatialRef()
# from Geometry
feature = layer.GetNextFeature()
geom = feature.GetGeometryRef()
spatialRef = geom.GetSpatialReference()
print spatialRef.ExportToWkt()

We choose to export to WKT here but doing dir(spatialRef) exposes various methods to export to different projection notations. It may help you to compare with other shapefiles.

See an excerpt from the output of dir(spatialRef)

 ...
 'ExportToMICoordSys',
 'ExportToPCI',
 'ExportToPrettyWkt',
 'ExportToProj4',
 'ExportToUSGS',
 'ExportToWkt',
 'ExportToXML',
 ...

In all case, you should also be aware that comparing projection is hard: the notation may change slightly but can affect or not the result you see on your screen or depending on the software (typically ArcGIS vs rest of the world)

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.