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 have 3 shape files with undefined projection, but I know which one it was used (Lisboa_Hayford_Gauss_IPCC). I need to 1st Define Projection. It works. I need to transform it from 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_1994 CS and it feets perfectly to the World Imagery provided by ArcGIS online or OSM basemap.

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

share|improve this question
    
1m is probably pushing the resolution and accuracy of the basemaps involved. Much of OSM is digitized from aerials with 1m or larger resolution. The World Imagery basemap is 0.3m resolution at best, 1m in most cases, and is built for a maximum scale of 1:70,000. On top of that, these basemaps are in Web Mercator Auxiliary Sphere. I don't think you should be concerned very much about a 1-2m deviation from the basemaps given all of these factors. –  blord-castillo Aug 12 '13 at 16:50

1 Answer 1

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.

share|improve this answer
    
I've tried the last transf. you suggested but the links now are more close to the side of the road instead of closer to the axis of the road (by looking to the OSM basemap and World Imagery Basemap). I don't know if this is better for me... do you have any suggestion? –  neverMind Aug 12 '13 at 16:41

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.