coordinate conversion, for excel/csv files

[JianLiu]

I have an excel file (or in fact converted to csv for easy use. Arcpy seems to be buggy about excel files...) with some lat/long coordinates that I would need to add corresponding projected coordinates. The way I can think of doing this is to use the MakeXYEvent to create an event layer, then reproject the layer and add the projected X and Y.

Implementation:

import os, arcpy, sys, traceback
from xlrd import open_workbook
import xlwt


arcpy.overwriteOutputs = True

projectionsPath = r"C:\Test\_UtilityScripts\Coordinate Systems_CopiedFromArcMap\Projected Coordinate Systems\State Plane\NAD 1983 (Meters)"
csvfile = r"C:\Test\6_USA_FIPs\excel\latlong_seperated\IL.csv"
lat_g = "lat_g"
long_g = "long_g"
outlayer = "xyEventlayer"
#savedlayer = r"C:\Test\6_USA_FIPs\python\temp\xyOutput.lyr"
#outputSHP = r"C:\Test\6_USA_FIPs\python\temp\xyOutput.shp"


if arcpy.Exists(outlayer):
    arcpy.Delete_management(outlayer)

srGCS = os.path.join(r"C:\Test\_UtilityScripts\Coordinate Systems_CopiedFromArcMap\Geographic Coordinate Systems\North America", "NAD 1983.prj")
srPRJ = os.path.join(projectionsPath, "NAD 1983 StatePlane Illinois East FIPS 1201 (Meters).prj")
arcpy.MakeXYEventLayer_management(csvfile, lat_g, long_g, outlayer, srGCS)

rows = arcpy.SearchCursor(outlayer, "", srPRJ)
row = rows.next()
feat = row.Shape
pnt = feat.getPart()
print pnt.X

I am using SearchCursor above with on the fly projection for coordinate conversion. But I am always getting coordinates of (nan, nan)??? I've tested the syntax with other shapefiles, and it should work... Anything special with this xyevent layer?  Please help, thanks people!

[MattEiben]

I'm a little unclear on what you're trying to do once you get the file into State Plane, but the process to getting it there could go something like this.  Note that here, I'm explicitly converting from WGS84 to State Plane:

arcpy.overwriteOutputs = True

csvfile = r"C:\Users\meiben\Desktop\test\fires.csv"
lat_g = "lat_g"
long_g = "long_g"
outlayer = "xyEventlayer"

arcpy.MakeXYEventLayer_management(csvfile, long_g, lat_g, outlayer,arcpy.SpatialReference(4326))
arcpy.CopyFeatures_management(outlayer,"templayer")
arcpy.Delete_management(outlayer)

arcpy.Project_management("templayer","coord_conv",arcpy.SpatialReference(26771))
arcpy.Delete_management("templayer")

A couple changes that I made here was using the Spatial Reference method instead of pointing to a file on the system.

outlayer,arcpy.SpatialReference(4326) = GCS WGS 84
arcpy.SpatialReference(26771) = PCS NAD_1927_StatePlane_Illinois_East_FIPS_1201

Check out the following URL for a list of all the codes:
http://resources.arcgis.com/en/help/main/10.1/018z/pdf/projected_coordinate_systems.pdf

At the end of this script, you should have a feature layer to do whatever you like with.

Hope this helps!

Matt

[JianLiu]

I'm a little unclear on what you're trying to do once you get the file into State Plane, but the process to getting it there could go something like this.  Note that here, I'm explicitly converting from WGS84 to State Plane:

arcpy.overwriteOutputs = True

csvfile = r"C:\Users\meiben\Desktop\test\fires.csv"
lat_g = "lat_g"
long_g = "long_g"
outlayer = "xyEventlayer"

arcpy.MakeXYEventLayer_management(csvfile, long_g, lat_g, outlayer,arcpy.SpatialReference(4326))
arcpy.CopyFeatures_management(outlayer,"templayer")
arcpy.Delete_management(outlayer)

arcpy.Project_management("templayer","coord_conv",arcpy.SpatialReference(26771))
arcpy.Delete_management("templayer")

A couple changes that I made here was using the Spatial Reference method instead of pointing to a file on the system.

outlayer,arcpy.SpatialReference(4326) = GCS WGS 84
arcpy.SpatialReference(26771) = PCS NAD_1927_StatePlane_Illinois_East_FIPS_1201

Check out the following URL for a list of all the codes:
http://resources.arcgis.com/en/help/main/10.1/018z/pdf/projected_coordinate_systems.pdf

At the end of this script, you should have a feature layer to do whatever you like with.

Hope this helps!

Matt

Thanks Matt for the reply! I wanted to use SearchCursor to do the projection on the fly (http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//000v00000039000000) and the code after the reprojection is just to print out the new coordinates for testing purposes, and that's where I got (nan,nan).

I wanted to use SearchCursor because it doesn't write anything to the harddrive, but the Project_management function would work too. I had just found out why my code failed -- it's not particularly due to the XYEventLayer, which I thought might give a problem to the on-the-fly projection. It is in fact due to bad coordinates in my file. The projection failed without giving an error.

As to the spatial reference input parameter, I think both the actual object and the string path would work. But it was also a factor that I had to filter out to pin down the problem. -- Sometimes arcpy is very confusing as to when it takes a string path or an actual object as input. There is no clear documentation at times, and it takes a long time to figure things out. Painful.

Thanks for the help, much appreciated!

[MattEiben]

Glad you figured out the problem!

I agree that it can be confusing sometimes as to what needs a string vs an object.  Though I've found that often they're interchangable in the arcpy.SpatialReference method.  For example, something like this would be totally valid:

WGS84 = """GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
        SPHEROID["WGS_1984",6378137.0,298.257223563]
    ],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]
]"""

arcpy.DefineProjection_management("featureClass",WGS84)

Here I literally opened up a .prj file and pasted the (formatted for readability) contents into my file. 

Also, I'm not sure which version of ArcMap you're running, but if you have 10.1+, you can use the new Data Access Cursors to get your coordinates, via the token shortcuts.  Something like this:

coords = [row[0] for row in arcpy.da.SearchCursor("FeatureClass",["SHAPE@XY"])]
>>> coords
>>> [(-117.54547981299241, 34.17563985328801)]

Though note that this method will only get the centroid, so if this was a polygon, you'd only get one point returned.  If you wanted to get all vertices from the geometry:

geom = [row[0] for row in arcpy.da.SearchCursor("FeatureClass",["SHAPE@"])][0]

partCount = geom.partCount

for n in range(partCount):
    part = geom.getPart(n)
    for point in part:
        print point.X, point.Y

Of course, the old school SearchCursor function still works great!  The new data access cursor methods are just a bit faster with larger datasets.

Matt