Projecting shapefile with transformation using OGR in python

Question

I want to change coordinate system of a shapefile from British National Grid to WGS84 (in order to be able to convert it to KML later on). To get the best results of projecting I would like to use OSGB_1936_To_WGS_1984_NGA_7PAR transformation (or later on define different transformation depending on location in UK).

The only description of this transformation that I know comes from ArcMap and looks like this:

Coordinate Frame - dx=446,000000 dy=-99,000000 dz=54,000000 rx=-0,9450000 ry=-0,261000 rz=0,435000 s=-20,892700

.. and I have no clue how to use it with osr.CoordinateTransformation(BNG, WGS84)

Any help highly appreciated!!

Jarek

Answer 1

Since OSR is based on Proj4 you are at the mercy of whatever it can accomplish. Have a look at this related post and the proj4 FAQ on towgs84. It appears from this that proj4 is handling the datum shift, etc. I would hope this would pass through into OGR/OSR, but can't comment specifically.

Assuming the other post is correct you could project your OGR features like so:

import ogr,osr
#... open an OGR feature layer

srcSR = osr.SpatialReference()
# Fixup the proj4 string as necessary for your projection
# or alternatively, find the appropriate EPSG code
srcSR.ImportFromProj4("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +units=m +no_defs")
destSR = osr.SpatialReference()
destSR.ImportFromEPSG(4326)  #wgs84

srTrans = osr.CoordinateTransformation(srcSR,destSR)

for ftr in ogrLyr:
    #Use the geometry inplace or clone it to decouple it from the layer
    newGeom = ftr.GetGeometryRef() #.Clone()

    #note, to get the actual coordinates you have to go one level deeper
    g = newGeom.GetGeometryRef(0)
    print(g.GetPoint(0))  #first point before transformation

    #apply the transformation
    newGeom.Transform(srTrans)

    print(g.GetPoint(0)) #same feature point after transformation