Calculating coordinates in different coordinate system programmatically

[TomBuckler]

I have a Python script that adds two fields and populates them with the X coordinate and Y coordinate of the centroid of the features. The feature class is in one projection, but I would like the centroid coordinates to be displayed in Lat/Lon. Is there a way to specify the coordinate system when you write something like !Shape.centroid!

I know when you right click on a field in the attribute table and select "Calculate Geometry" that you have the option to use the Data Frame's coordinate system. Is there a way to do this programmatically?

I've tried this function call and definition:

Lat(!Shape!)

def Lat(shape):
  arcpy.env.cartographicCoordinateSystem = "Coordinate Systems\Geographic Coordinate\Systems\World\WGS 1984.prj"
  lat = shape.centroid.x
  return lat

but it returns the x coordinate in the the same projection of the feature class rather than longitude.

[ChrisSnyder]

Yes, you just need to use an update cursor (and specify the spatial reference parameter):

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Setting_a_cursor_s_spatial_reference/0...

Then the "Shape" field's coordinate values will be reported in whatever coordinate system you specified...

[TomBuckler]

You sir are a gentlemen and a scholar. That's exactly what I was looking for.

If anyone is interested, here's my updated code for a script tool with a multivalue parameter.

import arcpy

latLonRef = "Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"

featureClasses = arcpy.GetParameterAsText(0)

featureClassesList = featureClasses.split(";")

for featureClass in featureClassesList:
 arcpy.AddMessage("Calculating XY coordinates for: " + featureClass)
 arcpy.AddField_management(featureClass, "LAT", "DOUBLE")
 arcpy.AddField_management(featureClass, "LON", "DOUBLE")
 rows = arcpy.UpdateCursor(featureClass, "", latLonRef)
 for row in rows:
  feat = row.shape
  coord = feat.getPart()
  lon = coord.X
  lat = coord.Y
  row.LAT = lat
  row.LON = lon
  rows.updateRow(row)
  #arcpy.AddMessage(str(lat) + ", " + str(lon))

[DaleHoneycutt]

At 10.0, there is a new tool, Convert Coordinate Notation, that does some interesting coordinate transformations.  It (probably) won't directly solve your exact problem, but I wanted to mention it for anyone stumbling upon this tread.

[danbecker]

You sir are a gentlemen and a scholar. That's exactly what I was looking for.

If anyone is interested, here's my updated code for a script tool with a multivalue parameter.

import arcpy

latLonRef = "Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"

featureClasses = arcpy.GetParameterAsText(0)

featureClassesList = featureClasses.split(";")

for featureClass in featureClassesList:
 arcpy.AddMessage("Calculating XY coordinates for: " + featureClass)
 arcpy.AddField_management(featureClass, "LAT", "DOUBLE")
 arcpy.AddField_management(featureClass, "LON", "DOUBLE")
 rows = arcpy.UpdateCursor(featureClass, "", latLonRef)
 for row in rows:
  feat = row.shape
  coord = feat.getPart()
  lon = coord.X
  lat = coord.Y
  row.LAT = lat
  row.LON = lon
  rows.updateRow(row)
  #arcpy.AddMessage(str(lat) + ", " + str(lon))

this is great; how would one modify to calc. area?

thanks!

[BerneJackson]

Couldn't get your code to work for polygon centroids.  Here's what I had to do to make it work:

import arcpy
latLonRef = "Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"
featureClasses = arcpy.GetParameterAsText(0)
featureClassesList = featureClasses.split(";")
for featureClass in featureClassesList:
arcpy.AddMessage("Calculating XY coordinates for: " + featureClass)
arcpy.AddField_management(featureClass, "LAT", "DOUBLE")
arcpy.AddField_management(featureClass, "LON", "DOUBLE")
rows = arcpy.UpdateCursor(featureClass, "", latLonRef)
for row in rows:
  feat = row.getValue("shape")
  cent = feat.centroid
                # To get the polygon area: cent = feat.area
  row.LAT = cent.Y
  row.LON = cent.X
  rows.updateRow(row)
  #arcpy.AddMessage(str(lat) + ", " + str(lon))

[BobNoyes]

This works for what I wanted except it locks up ArcMap 10.2.2 when I run it. When I check the attribute table the columns are created and populated.  Any ideas why this would "freeze up" ArcMap and not continue the other request in the script?

Any suggestions would be appreciated!

Addendum: Okay, I wasn't patient enough. It takes almost 30 minutes for the lat long function to run. It will then complete the python script. Not sure why it takes so long. I am wondering if it has something to do with it being in a Sql GDB? After testing the lat long it looks like they are in the GDB projection (state plane) not WGS84, which is what I wanted. I am going to try Plan B and set up a staging file to add the fields and calculations in a shapefile then copy to the GDB.

[danbecker]

Hi everyone,

I have a gdb (v9.3) with several feat. datasets. All the feat. datasets are unprojected, coor. system set to WGS1984. I can calc. coord_x and coord_y fields of the point feat. classes fine using the below code, all is well.

Once the script hits a polygon feat. class, I would like to calc. the area_size field to area - acres. At present, the code doesn't calculate anything, even with the polyref pointing to a projected, .prj file.

At the bottom, in the 'polygon' section, you'll see a commented out section:

#expression = "!shape.area@acres!"                                    
 #gp.CalculateField_management(fc, fn, expression, "PYTHON_9.3")

When this is uncommented, (and the updatecursor method commented out) it calculates area just fine, but of course, since the feat. class is unprojected, it calc.s area_size to 1.4857433E-06, which obiv. is in WGS84 DD.

Any thoughts?

thanks!

# Import system modules
import arcgisscripting, os

gp = arcgisscripting.create(9.3)

# Load required toolboxes...
#gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
gp.AddToolbox("C:/Program Files/ArcGIS/Desktop10.0/ArcToolbox/Toolboxes/Data Management Tools.tbx")

#gp.workspace = r"G:\\63D RRC GIS Data\\63D_RSC_GDB"
gp.workspace = r"C:\\temp\\python"
#gp.workspace = r"G:\\63D RRC GIS Data"

pathwrksp = gp.workspace
# print pathwrksp
PointRef = r"C:\\Program Files\\ArcGIS\\Desktop10.0\\Coordinate Systems\\Geographic Coordinate Systems\\World\\WGS 1984.prj"
PolyRef = r"C:\\Program Files\\ArcGIS\\Desktop10.0\\Coordinate Systems\\Projected Coordinate Systems\\Continental\\North America\\USA Contiguous Albers Equal Area Conic USGS.prj"


path2wrksps = gp.ListWorkspaces("OK002*")

for path2wrksp in path2wrksps:
    #print "... " + path2wrksp
    gp.workspace = path2wrksp
    facil_id = str(path2wrksp[34:39])
    instln_id = "06510"
    #print facil_id
    wrksps3 = gp.ListWorkspaces("*", "FileGDB")
    for wrksp3 in wrksps3:
        gp.workspace = wrksp3
        print "...... " + gp.workspace
        datasets = gp.ListDatasets("cad*")        
        for dataset in datasets:
            gp.workspace = wrksp3 + "\\" + dataset
            #print dataset
            FcList = gp.ListFeatureClasses("*","Point")
            for fc in FcList:                
                gp.workspace = dataset + "\\" + fc
                if fc <> "soil_map_unit_area":
                    if fc <> "wetland_area":
                        if fc <> "flood_zone_area":
                            print fc
                            #if gp.exists("vehicle_parking_area"):
                            rows = gp.UpdateCursor(fc, "", PointRef)
                            for row in rows:
                                allFields = gp.listfields(fc)
                                print fc
                                for field in allFields:
                                    fn = field.name
                                    if fn == "coord_x":
                                        feat = row.shape
                                        coord = feat.getPart()
                                        lon = coord.X
                                        row.coord_x = lon
                                        rows.updateRow(row)                                        
                                        print fc, fn + " calculated to long"
                                    elif fn == "coord_y":
                                        feat = row.shape
                                        coord = feat.getPart()
                                        lat = coord.Y
                                        row.coord_y = lat
                                        rows.updateRow(row)
                                        print fc, fn + " calculated to lat"
            FcList = gp.ListFeatureClasses("*","Polygon")
            for fc in FcList:                
                gp.workspace = dataset + "\\" + fc
                if fc <> "soil_map_unit_area":
                    if fc <> "wetland_area":
                        if fc <> "flood_zone_area":
                            print fc
                            rows = gp.UpdateCursor(fc, "", PolyRef)
                            for row in rows:
                                allFields = gp.listfields(fc)
                                print fc
                                for field in allFields:
                                    fn = field.name
                                    if fn == "area_size":
                                        feat = row.getValue("SHAPE")
                                        cent = feat.area
                                        rows.updateRow(row)
                                        #expression = "!shape.area@acres!"                                    
                                        #gp.CalculateField_management(fc, fn, expression, "PYTHON_9.3") 
                                        print fc, fn + " calculated to area"

[RobinSmith]

Our geodatabases are set to State Plane projection and the x, y, z fields display coordinates in that projection.  I've successfully used Calculate Geometry to populate the State Plane X, Y fields from points created in ArcMap. 
For reporting purposes, I also need to have Lat/Lon fields populated as well for each of those same features (point).  How can I use Calculate Field, or Calculate X, Y, or Add X, Y to take the State Plane coordinates in those fields to calculate and populate the Lat/Lon fields?  I'm not a programmer, so step by step instructions would be very much appreciated.  Thanks.

[danbecker]

depending on how many gdb's you have, this can be done in ArcMap...manually.

1. Start with a new ArcMap document, then double click on 'layers' in the table of contents
2. in data frame properties, go to coordinate system tab/ predefined/ geographic coor system/world/wgs84
3. hit "ok, then "yes" on the warning
4. add your feat classes
5. create 2 new fields, named 'lat' and 'long'
6. right click on the field and choose 'calc geometry'
7. select 'y coordinate of a point' then select 'use coordinate system of the dataframe' units 'decimal degrees'
8. you now have a column of latitude coordinates.
9. do the same for 'long' field.

remember, x measures distance E/W of Prime Meridian, y measures N/S of equator