Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question

1 Answer 1

up vote 5 down vote accepted

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
share|improve this answer
    
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 2: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 2:57

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.