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

score 6 · Accepted Answer · edited May 23 '17 at 12:39

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

Thanks tharen, that helps a lot. It seems like the transformation that I have from ArcMap reflects seven parameters notation delta_x, delta_y, delta_z, Rx - rotation X, Ry - rotation Y, Rz - rotation Z, M_BF - Scaling, so my definition of British National Grid with NGA 7PAR transformation in proj4 would look like this: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +towgs84=446.0,-99.0, 54.0,-0,9450000,-0.261000,0.435000,-20,892700 +units=m +no_defs — jareks, Jan 28 '12 at 02:11
I have seen discussions on the web comparing ArcGIS projections to proj4, they can be slightly different. It looks like you are on the right track though. If you can get a hold of some reference data that has been projected to both coordinate systems you could double check against that. — tharen, Jan 28 '12 at 02:57

Projecting shapefile with transformation using OGR in python

1 Answers1

Linked