1

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?

3
  • Its likely there's more than one AUTHORITY attribute within the proj variable. Try inspecting the value of the proj variable as WKT to see what's going on. Compare how it looks for one of the problem new images versus an older image for which the code worked.
    – M Bain
    Aug 18 '21 at 11:55
  • @MBain I have added the getProjection output. Am I to believe getProjection isn't perfect and might not correctly create the projection?
    – g0atm1lk
    Aug 18 '21 at 12:21
  • Not exactly sure g0atm1lk, I've just been grappling with a not dissimilar situation - trying to determine the EPSG code for a number of vector layers, and found the contents of SpatialReference varied greatly between files. GetAttrValue('AUTHORITY',1) wasn't always the AUTHORITY value I was looking for. See user2856's answer - look's informative.
    – M Bain
    Aug 18 '21 at 20:49
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

enter image description here

1
  • Thank you. It works a charm.
    – g0atm1lk
    Aug 19 '21 at 10:11

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.