in python, srs.ExportToProj4() returns an empty string, why?

Question

I have a shapefile with the following .prj file:

PROJCS["NAD_1927_Albers",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],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",23.0],UNIT["Meter",1.0]]

When I run the following python code:

    srs = osr.SpatialReference()
    if ( srs.ImportFromWkt( prjText ) != 0 ):
        return "error"

    srsProj4 = srs.ExportToProj4()
    sourceProjection = pyproj.Proj( srsProj4, preserve_units=True )

the variable srsProj4, as returned from srs.ExportToProj4(), is an empty string.

When I look here: http://www.gdal.org/ogr/classOGRSpatialReference.html

it says that ExportToProj4() returns an empty string in the case the srs is in a local (non-geo-referenced) coordinate system. But as far as I understand, the text in the .prj file (above) is a valid, regular geo projection.

This code works with other .prj files I've used. Can you tell what's special about this one that's making it fail?

Thanks.

Answer 1

It seems that osr has difficulties with PROJECTION["Albers"]. This could be a bug, and it is odd that gdalsrsinfo interprets, and corrects this to PROJECTION["Albers_Conic_Equal_Area"].

Here's the hackish fix:

from osgeo import osr
osr.UseExceptions()

prjText = 'PROJCS["NAD_1927_Albers",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],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",23.0],UNIT["Meter",1.0]]'

prjText = prjText.replace('"Albers"', '"Albers_Conic_Equal_Area"')
srs.ImportFromWkt(prjText)  # this will now raise RuntimeError for corrupt data
print(srs.ExportToProj4())  # another RuntimeError will be raised on failure

now shows:

+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs

Answer 2

gdalsrsinfo recognizes the definition, returning:

PROJ.4 : '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +
datum=NAD27 +units=m +no_defs '

OGC WKT :
PROJCS["NAD_1927_Albers",
    GEOGCS["GCS_North_American_1927",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke_1866",6378206.4,294.9786982]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["longitude_of_center",-96.0],
    PARAMETER["Standard_Parallel_1",29.5],
    PARAMETER["Standard_Parallel_2",45.5],
    PARAMETER["latitude_of_center",23.0],
    UNIT["Meter",1.0]]

which should be equivalent to EPSG:5069 named NAD27 / Conus Albers

Maybe the name is not recognized correctly in Python.