topic Calculating coordinates in different coordinate system programmatically in Geoprocessing Questions

Calculating coordinates in different coordinate system programmatically

TomBuckler — Fri, 19 Nov 2010 16:33:55 GMT

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.

Re: Calculating coordinates in different coordinate system programmatically

ChrisSnyder — Fri, 19 Nov 2010 17:46:53 GMT

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/002z0000001w000000/

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

Re: Calculating coordinates in different coordinate system programmatically

TomBuckler — Sat, 11 Dec 2021 16:03:04 GMT

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))

Re: Calculating coordinates in different coordinate system programmatically

DaleHoneycutt — Fri, 19 Nov 2010 20:42:05 GMT

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.

Re: Calculating coordinates in different coordinate system programmatically

danbecker — Sat, 11 Dec 2021 16:03:06 GMT

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!

Re: Calculating coordinates in different coordinate system programmatically

BerneJackson — Tue, 29 Mar 2011 21:44:32 GMT

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))

Re: Calculating coordinates in different coordinate system programmatically

danbecker — Sat, 11 Dec 2021 16:03:09 GMT

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"

Re: Calculating coordinates in different coordinate system programmatically

RobinSmith — Wed, 27 Apr 2011 20:18:41 GMT

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.

Re: Calculating coordinates in different coordinate system programmatically

danbecker — Fri, 03 Jun 2011 19:53:16 GMT

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

Re: Calculating coordinates in different coordinate system programmatically

danbecker — Fri, 03 Jun 2011 20:00:02 GMT

not the most elegant solution, but to calc. geometries of an unprojected feat. class, I ended up:

projecting feat. class - output to temp.gdb
calculated geometry of temp feat. class
make feat. layer of temp feat. class
join temp feat. layer "OBJECTID" ----- with ----- original feat. class "OBJECTID"
calc field of orig. feat. class with the calculated geometries

Re: Calculating coordinates in different coordinate system programmatically

MichaelWalden — Wed, 18 Dec 2013 19:15:36 GMT

I am very new to Python and have a similar challenge.

I have a point feature class with about 8,500 records that is in WGS84. The feature class already has 6 fields (data type = double, format = Decimal Degrees) created to calculate the longitude and latitude in 3 different coordinate systems (WGS84_Longitude, WGS84_Latitude, NAD83_Longitude, NAD83_Latitude, NAD27_Longitude, and NAD27_Latitude). I am currently manually doing this in an ArcMap edit session, calculating the geometry using the dataframe's coordinate system, however I would like to automate the process and schedule it to run at a predetermined frequency.

I need to write a python script that will calculate the 6 fields using the correct coordinate system and a predetermined datum transformation. Can this be done?

Re: Calculating coordinates in different coordinate system programmatically

BobNoyes — Tue, 17 Apr 2018 15:57:02 GMT

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.