How to assign OGR geometry a spatial reference when writing to a new layer using the Python API?

Question

I'm using the Python OGR API to add new features to a PostGIS layer. However, despite setting the geometry spatial reference, when the geometry is written out the spatial reference is missing. This is a problem as the destination table has check constraint forcing geom_srid to be a particular type.

For example:

ogr.UseExceptions()

# Create a feature in an existing layer (where layer is a PostGIS table)
feature = ogr.Feature(layer.GetLayerDefn())

# Create spatial reference
out_srs = ogr.osr.SpatialReference()
out_srs.ImportFromEPSG(27700)

# Assign to geometry
geom.AssignSpatialReference(out_srs)

# Set feature geometry
feature.SetGeometry(geom)

# Create the feature
layer.CreateFeature(feature)
feature.Destroy()

Returns:

Error: new row for relation "table" violates check constraint 
"enforce_srid_geom"
Command: INSERT INTO "table" ("geom", "attr1") VALUES (GeomFromEWKT('SRID=-1;POINT (5000, 1000)'::TEXT), 'row1')

At this point it is clear that OGR is converting the geometry to an EWKT string without pulling in the spatial reference attribute.

Specific questions:

1) Does anyone know how I can force the geometry output from SetGeometry() to include an SRID when it is converted to EWKT?

2) Has anyone found this with the CPP bindings?

3) Is this a bug in OGR? From OGR source the output string is being created in ogrpgtablelayer.cpp where nSRSID is being set to -1 if no valid SRS input is found.

Many thanks,

Tom

Answer 1

Your approach seems correct and I am not sure where the problem is. As an alternate, you can try using ImportFromProj4() and supply a proj4 string for EPSG 27700 and see if that works.