How to read geographic coordinates when the Shapefile has a projected spatial reference

Question

I have a polygon Shapefile and want to read the polygon's vertices using OGR and Python. The reference system of the file is projected. However, for further processing I need the geographic coordinates of the vertices.

This is my code:

import ogr
path_to_infile = '/path/to/shapefile.shp' 
driver = ogr.GetDriverByName('ESRI Shapefile') 
infile = driver.Open(path_to_infile) 
layer = infile.GetLayer()
for area in layer: 
    area_shape = area.GetGeometryRef() 
    area_polygon = area_shape.GetGeometryRef(0) 
    no_of_polygon_vertices = area_polygon.GetPointCount() 
    for vertex in xrange(no_of_polygon_vertices): 
        lon, lat, z = area_polygon.GetPoint(vertex)

Coordinates (lat & lon variables) will be read according to the spatial reference. In my case, projected coordinates are displayed. Is there a way to directly extract geographic coordinates in OGR?

I know that .prj files of projected reference systems also store information on the corresponding geographic coordinate system. In arcpy for instance, geographic coordinates of a Shapefile with projected spatial reference can be accessed using the GCS property of arcpy.Describe(dataset).spatialReference:

import arcpy
# accessing projected coordinates
spatial_ref = arcpy.Describe(dataset).spatialReference 

# accessing geographic coordinates
spatial_ref.GCS

When spatial_ref.GCS is defined, I get geographic point coordinates for a polygon's vertices.

Is there a similar property for OGR available? Reprojecting is not my first choice, as I'm working in a custom reference system.

Answer 1

For you, what is the difference between geographic coordinates and projected coordinates ?

From ESRI GeoNet: Projected Coordinate System vs. Geographic Coordinate System:

• A Geographic Coordinate Systems is defined by a 3-D surface and is measured in latitude and longitude.
• A Projected Coordinate systems refers to data that is defined by a flat 2-D surface and can be measured in units of meters and feet.
• When displaying data that's using a geographic coordinate system, ArcMap uses a 'pseudo-Plate Carree' projection. Basically, it just treat the coordinate values as if they're linear and just display the data.

Therefore if your coordinates are lat & lon variables, you use a Geographic Coordinate System (sensu ESRI)

1) With osgeo.ogr and a shapefile with WGS84 projection (Geographic Coordinate System), the crs is found by:

from osgeo import ogr
infile = ogr.Open("aWGS84.shp")
layer = infile.GetLayer()
# crs
spatialRef = layer.GetSpatialRef()
spatialRef.ExportToWkt()
'GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]'

You can also use spatialRef.ExportToPrettyWkt(),spatialRef.ExportToProj4(),spatialRef.ExportToUSGS() and many other formats.

First feature of the shapefile

feature = layer.GetFeature(0)
feature.ExportToJson()
'{"geometry": {"type": "Polygon", "coordinates": [[[5.826429692987699, 50.63093772365397], [5.830982302422351, 50.62928355089628], [5.826407932565759, 50.6268451632606], [5.82334832264974, 50.62941225500874], [5.826429692987699, 50.63093772365397]]]}, "type": "Feature", "properties": {"id": null}, "id": 0}'

And the point coordinates of the Polygon (WGS84 in degrees)

import json
json.loads(feature.ExportToJson())['geometry']['coordinates']
[[[5.826429692987699, 50.63093772365397], [5.830982302422351, 50.62928355089628], [5.826407932565759, 50.6268451632606], [5.82334832264974, 50.62941225500874], [5.826429692987699, 50.63093772365397]]]

2) With the same shapefile in a projected coordinate system (in meters)

infile = ogr.Open("aproj.shp")
layer = infile.GetLayer()
feature = layer.GetFeature(0)
json.loads(feature.ExportToJson())['geometry']['coordinates']
[[[253122.92278611325, 147711.88454888575], [253448.57840302447, 147534.25421241205], [253130.3240500753, 147256.7068115687], [252908.28612947493, 147537.9548443528], [253122.92278611325, 147711.88454888575]]]

3) For the conversion between the two geometries, I need to reproject the shapefile (degrees to meters) with osgeo.osr or Pyproj (there are many examples in GIS SE).

4) Instead of using osgeo.ogr why don't use the more "Pythonic" Fiona (another Python wrapper of the OGR library)

import fiona
file = fiona.open("aWGS84.shp")
spatialRef = file.crs
spatialRef
{'init': u'epsg:4326'}

Fiona use Python dictionaries for all (GeoJson format), therefore the first feature is obtained by

file.next()
{'geometry': {'type': 'Polygon', 'coordinates': [[(5.826429692987699, 50.63093772365397), (5.830982302422351, 50.62928355089628), (5.826407932565759, 50.6268451632606), (5.82334832264974, 50.62941225500874), (5.826429692987699, 50.63093772365397)]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', None)])}

Answer 2

In order to get geographic coordinates from a Shapefile with projected spatial reference, the geometries need to be reprojected when working with OGR. I tried to avoid this at first, because I'm working in a custom spatial reference system and I'm not sure about the reprojection parameters. When I look at the reprojected data in a GIS, they seem to be accurate, though.

I added the following lines to the script:

sr = layer.GetSpatialRef() 
geogr_sr = sr.CloneGeogCS()
proj_to_geog = osr.CoordinateTransformation(sr, geogr_sr)
area_polygon.Transform(proj_to_geog)

The final code looks like this:

import ogr
path_to_infile = '/path/to/shapefile.shp' 
driver = ogr.GetDriverByName('ESRI Shapefile') 
infile = driver.Open(path_to_infile) 
layer = infile.GetLayer()
# get projected spatial reference
sr = layer.GetSpatialRef()
# get geographic spatial reference 
geogr_sr = sr.CloneGeogCS()
# define reprojection
proj_to_geog = osr.CoordinateTransformation(sr, geogr_sr)
for area in layer: 
    area_shape = area.GetGeometryRef() 
    area_polygon = area_shape.GetGeometryRef(0) 
    # transform coordinates
    area_polygon.Transform(proj_to_geog)
    no_of_polygon_vertices = area_polygon.GetPointCount() 
    for vertex in xrange(no_of_polygon_vertices): 
        lon, lat, z = area_polygon.GetPoint(vertex)

The differences in speed from the reprojection are neglectable in my opinion (< 0.01 sec on average with ~35000 vertices on an off-the-shelf laptop).

This post was very helpful to me.