[GIS] Geographic Transformations using ArcPy vs. Coordinate Systems used by ArcGis online basemaps

arcpycoordinate system

I have three shapefiles with undefined projection, but I know which one was used (Lisboa_Hayford_Gauss_IPCC). I need to first Define Projection. It works. I then need to transform it to WGS_1984, so I've used:

import sys, arcpy
#globals set env.workspace and other used paths
import globals as g


def main():
    # set local variables
    inDataset = ["rede_Coimbra_node.shp", "rede_Coimbra_link.shp", "rede_Coimbra_zone.shp"]

    prjFile = g.coveragePath + "Lisboa Hayford Gauss IPCC.prj"
    prjFileDatum73 = g.coveragePath + "Datum 73 Hayford Gauss IPCC.prj"
    prjFileWGS = g.coveragePath + "WGS 1984.prj"
    # Geographic transformation - 
    transformation1 = "Datum_Lisboa_Hayford_To_Datum_73_1"
    transformation2 = "Datum_73_To_WGS_1984_4"

    sr1 = arcpy.SpatialReference(prjFile)
    print sr1.name
    #"spatialRef.name = Lisboa_Hayford_Gauss_IPCC"
    #sr2= arcpy.SpatialReference(prjFileDatum73)
    #print sr2.name
    #print (arcpy.ListTransformations(sr1, sr2))
    #spatialRef.name = "Datum_73_Hayford_Gauss_IPCC"
    sr3 = arcpy.SpatialReference(prjFileWGS)
    print (arcpy.ListTransformations(sr1, sr3))
    #spatialRef.name = "GCS_WGS_1984"

    try:
        print "..."
        for fc in inDataset:
            fcFullPath = g.coveragePath+fc
            #arcpy.DefineProjection_management(fcFullPath, prjFile)
            outfc= fc.strip(".shp")
            outfc = arcpy.env.workspace+"/"+outfc+"_wgs1984"
            print outfc

            arcpy.Project_management (fcFullPath, outfc,prjFileWGS, u'Datum_Lisboa_Hayford_To_Datum_73_1 + Datum_73_To_WGS_1984_4',prjFile)          
            #arcpy.Project_management (fcFullPath, outfc,prjFileDatum73, transformation1,prjFile)
            #arcpy.Project_management (outfc, outfc+"_wgs84",prjFileWGS, transformation2,prjFileDatum73)

    except Exception, e:
        # If an error occurred, print line number and error message
        tb = sys.exc_info()[2]
        print "-->Line %i" % tb.tb_lineno
        print e.message

if __name__ == '__main__':
    main()

I also have a layer from OSM, which I also projected to WGS_1984 GCS and it fits perfectly to the World Imagery provided by ArcGIS Online or OSM basemap.

However, my projected three shapefiles (now converted to feature classes) don't seem to fit so well. They have a deviation from 1-2m to the raster basemap. Did I do something wrong or is it because the online basemaps are in other CRSs and I can't change that?

Best Answer

If you have ArcGIS 10.1, I would change the second transformation to "Datum_73_To_WGS_1984_2009_7par". It uses the same parameters as EPSG:5037, Datum 73 to ETRS89 (5). That transformation is listed has having an accuracy of 2m, but the accuracy should still be better than the one you're using. There's also a more accurate 7 parameter transformation for Datum Lisboa Hayford to Datum 73, "Datum_Lisboa_Hayford_To_Datum_73_2" which has a listed accuracy of 2m in comparison to 5m for "Datum_Lisboa_Hayford_To_Datum_73_1".

There's also a 7 parameter transformation directly from Datum Lisboa Hayford to WGS84, "Datum_Lisboa_Hayford_To_WGS_1984_2009_7par" with an accuracy around 1.5m to 2m.

Related Question