GDAL extracted EPSG not compatible with Fiona

Question

I am creating Shapefiles from features extracted from JP2 images in python. My current method of creating Shapefiles by extracting the EPSG from the source image with GDAL and using Fiona is as below and works for most images.

# getEPSG from image
ds = gdal.Open("jp2_image.jp2", gdal.GA_ReadOnly)
proj = osr.SpatialReference(wkt=ds.GetProjection())
epsg = proj.GetAttrValue('AUTHORITY',1)

# create shape file
with fiona.open(shp_file_output, 'w', 'ESRI Shapefile', schema=schema, crs=from_epsg(espg)) as shp_file:
        shp_file.write({
            'geometry' : {'type':'MultiPolygon',
            'coordinates': [multi_polygon]},
            'properties': {
                'tag': 'Image features'
            }
        })

However, I received some new images and Fiona throws an error with the extracted EPSG.

Invalid input to create CRS: {'init': 'epsg:6326', 'no_defs': True}

What is interesting is I can see the correct EPSG in the image metadata, the EPSG being 32756. When I manually enter 32756 the shapefiles are correctly generated.

with fiona.open(shp_file_output, 'w', 'ESRI Shapefile', schema=schema, crs=from_epsg(32756)) as shp_file:

Here is the get projection output for the file that works:

'PROJCS["GDA2020 / MGA zone 56",GEOGCS["GDA2020",DATUM["Geocentric_Datum_of_Australia_2020",SPHEROID["GRS 1980",6378137,298.257222100887,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1168"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","7844"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","7856"]]'

and for file that doesn't, no mention of 32756:

'PROJCS["WGS_84_UTM_zone_56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'

How do I get GDAL to return the correct EPSG needed to locate the Shapefiles correctly?

Answer 1

The authority code you're extracting isn't always the code for the whole coordinate system. In this case, GDAL is giving you the code for the WGS84 datum EPSG:6326 (not the WGS84 Geodetic CRS which is EPSG:4326).

Simplest way of doing this is to pass the wkt to fiona.open and let fiona generate a CRS.

Basic example:

import fiona

# wkt = ds.GetProjection()  # commented out just to demo, hardcoded wkt below is the same as yours.
wkt = 'PROJCS["WGS_84_UTM_zone_56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
schema = {'geometry': 'Point', 'properties': {'id': 'int'}}

with fiona.open('/tmp/test.shp', 'w', 'ESRI Shapefile', schema, crs=wkt) as layer:
    print(layer.crs)  # I'm creating an empty file here just to  demo...

{'init': 'epsg:32756'}

And if you use rasterio instead of GDAL, you can just pass the raster dataset crs property to fiona:

import fiona
import rasterio

with rasterio.open('your.jp2') as src:
    crs = src.crs

schema = {'geometry': 'Point', 'properties': {'id': 'int'}}

with fiona.open('/tmp/test.shp', 'w', 'ESRI Shapefile', schema, crs=crs) as layer:
    pass