xangis

Trouble With Snap_edit Data

Discussion created by xangis on May 8, 2014
I have a Python script that takes a map file containing a hydrography layer, loads a CSV file into a layer containing a set of points, and then snaps those points to the hydrography layer. After doing that, it exports the data back to a CSV file. (The CSV file will be used in another application).

When I open the .mxd file in ArcMap 10.2, everything looks good, the points are snapped as expected.

The trouble is, when I export the data, it's in the pre-snapped form without the Snap transformation applied. I can confirm visually this by adding the output CSV data to the map in ArcMap. Is there some other method I should be using to export the snapped data layer as a CSV instead of arcpy.SearchCursor?

A second issue I'm having is that the Project_management command doesn't appear to actually do anything to the coordinate data as far as I can tell. Is there a command I can use to translate coordinates back and forth between WGS1984 and NAD1983? I'm OK with translating them point by point if need be.

Here's my code:
def SnapSamplesToHydrography():
    """
    Snaps samples in a CSV to a hydrography layer.
    """
    print 'Starting SnapSamplesToHydrography'
    # Import system modules
    import arcpy
    import arcpy.mapping as mapping
    # Set local variables
    input_csv = "C:/Sitka/samplesites.csv"
    with open(input_csv, 'r') as csv_in:
        rows = csv.reader(csv_in, quoting=csv.QUOTE_NONE)
        for row in rows:
            print 'Row: ' + str(row)
        csv_in.close()
    arcpy.env.workspace = "C:/Sitka/shp/"
    mapfile = "C:/Sitka/Hydrography.mxd"
    working_mapfile = "C:/Sitka/Hydrography_Active.mxd"
    reproj_layer = "C:/Sitka/Hydrography_Active_Reproj.shp"
    copied_samples = "C:/Sitka/SampleSitesCopy.shp"
    mxd = mapping.MapDocument(mapfile)
    mxd.saveACopy(working_mapfile)
    mxd = mapping.MapDocument(working_mapfile)
    print "Map Layers: "
    layers = mapping.ListLayers(mxd)
    print layers
    print "Map Data Frames: "
    frames = mapping.ListDataFrames(mxd)
    #print frames
    for frame in frames:
        print("Data Frame: " + frame.name)
    for borken in mapping.ListBrokenDataSources(mxd):
        print("Borken: " + borken.name)
    # Creating a spatial reference object
    sr = arcpy.SpatialReference()
    sr.factoryCode = 3857
    sr.create()
    othersr = arcpy.SpatialReference()
    othersr.factoryCode = 4269
    othersr.create()
    arcpy.MakeXYEventLayer_management(input_csv, "Field3", "Field4", "samplesites", sr)
    print("MakeXYEventLayer result: " + arcpy.GetMessages())
    print("SaveToLayerFile result: " + arcpy.GetMessages())
    result = arcpy.CopyFeatures_management("samplesites", copied_samples)
    #arcpy.MakeFeatureLayer_management(copied_samples, samples_layer)
    layer = mapping.Layer(copied_samples)
    mapping.AddLayer(frames[0], layer)
    mxd.save()
    print("CopyFeatures_management result: " + str(result))
    print("CopyFeatures_management result status: " + str(result.status))
    print("CopyFeatures_management getOutput: " + result.getOutput(0))
    print("Map Layers: ")
    layers = mapping.ListLayers(mxd)
    print(layers)
    result = arcpy.Snap_edit(layers[0], [[layers[1], "EDGE", "500 Meters"]])
    print("Snap result: " + arcpy.GetMessages())
    print("Return Value: " + str(result))
    print('Result status: ' + str(result.status))
    print(result.getOutput(0))
    result = arcpy.Project_management(layers[0], reproj_layer, othersr, transform_method="WGS_1984_(ITRF00)_To_NAD_1983")
    print(result)
    print('Result status: ' + str(result.status))
    print(result.getOutput(0))
    print "Map Layers: "
    layers = mapping.ListLayers(mxd)
    print layers
    desc = arcpy.Describe(working_mapfile)
    print("Describe: " + str(desc))
    print("Describe dir: " + str(dir(desc)))
    with open('C:/Sitka/out.csv', 'w') as csv_out:
        rows = arcpy.SearchCursor(layers[0])
        writer = csv.writer(csv_out, quoting=csv.QUOTE_NONE)
        for row in rows:
            print("Full Row: " + str(row.Field1) + ", " + str(row.Field2) + ", " + str(row.Field3) + ", " +
                  str(row.Field4) + ", ")
            writer.writerow([row.Field1, row.Field2, row.Field3, row.Field4])
        csv_out.close()
    print arcpy.GetMessages()
    print 'Finished SnapSamplesToHydrography'


This is the text output for the command:

Starting SnapSamplesToHydrography
Row: ['2190247', '1053', '-12782895.341086963', '5527399.948853009']
Row: ['2190248', '1053', '-12771785.652828144', '5520218.6609102385']
...
Row: ['2190335', '1053', '-12771396.034772772', '5522220.941711186']
Row: ['2190336', '1053', '-12770919.586732037', '5517841.71775302']
Map Layers:
[<map layer u'NHDFlowline.lyr'>]
Map Data Frames:
Data Frame: Layers
Frame zero: <geoprocessing Data Frame object object at 0x19D972A0>
Tables:
MakeXYEventLayer result: Executing: MakeXYEventLayer C:/Sitka/samplesites.csv Field3 Field4 samplesites "PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]];-20037700 -30241100 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision" #
Start Time: Thu May 08 15:29:28 2014
Succeeded at Thu May 08 15:29:28 2014 (Elapsed Time: 0.01 seconds)
SaveToLayerFile result: Executing: MakeXYEventLayer C:/Sitka/samplesites.csv Field3 Field4 samplesites "PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]];-20037700 -30241100 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision" #
Start Time: Thu May 08 15:29:28 2014
Succeeded at Thu May 08 15:29:28 2014 (Elapsed Time: 0.01 seconds)
CopyFeatures_management result: C:\Sitka\SampleSitesCopy.shp
CopyFeatures_management result status: 4
CopyFeatures_management getOutput: C:\Sitka\SampleSitesCopy.shp
Map Layers:
[<map layer u'SampleSitesCopy'>, <map layer u'NHDFlowline.lyr'>]
Snap result: Executing: Snap GPL0 "GPL1 EDGE '500 Meters'"
Start Time: Thu May 08 15:29:28 2014
Succeeded at Thu May 08 15:29:41 2014 (Elapsed Time: 13.26 seconds)
Return Value: GPL0
Result status: 4
GPL0
C:\Sitka\Hydrography_Active_Reproj.shp
Result status: 4
C:\Sitka\Hydrography_Active_Reproj.shp
Map Layers:
[<map layer u'SampleSitesCopy'>, <map layer u'NHDFlowline.lyr'>]
Describe: <geoprocessing describe data object object at 0x19D90F38>
Describe dir: []
Full Row: 2190247, 1053, -12782895.3411, 5527399.94885,
Full Row: 2190248, 1053, -12771785.6528, 5520218.66091,
...
Full Row: 2190335, 1053, -12771396.0348, 5522220.94171,
Full Row: 2190336, 1053, -12770919.5867, 5517841.71775,
Executing: Project GPL0 C:\Sitka\Hydrography_Active_Reproj.shp GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]] WGS_1984_(ITRF00)_To_NAD_1983 PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]]
Start Time: Thu May 08 15:29:43 2014
Succeeded at Thu May 08 15:29:43 2014 (Elapsed Time: 0.14 seconds)
Finished SnapSamplesToHydrography


Please be gentle. I'm not new to Python, but I am very new to ArcGIS.

Outcomes