Projecting shapefile with transformation using OGR in python

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

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