Getting latitude and longitude of projected point using ArcPy? [closed]

Question

I have a point feature in a feature class that is being accessed by ArcPy. The point is projected but I need to find an efficient means to get the unprojected latitude and longitude for this point.

Is there a method other than reprojecting (unprojecting), getting a search cursor on the new feature class, finding the feature, then getting the lat/lon off the feature's shape?

Answer 1

Most of the other answers were all posted when ArcGIS 10.0 was the latest software. At ArcGIS 10.1 much new ArcPy functionality became available. This answer takes advantage of that new functionality. It will not be suitable for 10.0 but offers increased performance and functionality for 10.1 and later.

import arcpy

input_feature_class = 'C:\your_feature_class.shp'
wkid = 4326 # wkid code for wgs84
spatial_reference = arcpy.SpatialReference(wkid)

fields_to_work_with = ['SHAPE@']

with arcpy.da.SearchCursor(input_feature_class,
                           fields_to_work_with) as s_cur:
    for row in s_cur:
        point_in_wgs84 = row[0].projectAs(spatial_reference)
        print point_in_wgs84.firstPoint.X, point_in_wgs84.firstPoint.Y

This code snippet uses the wkid to create a spatial reference object rather than typing out a string representation, uses the more modern data access cursors, and projects the individual geometry objects using the projectAs() method.

Answer 2

The SearchCursor supports specifying a spatial reference- in this case, you'd want a Geographic Coordinate System, such as WGS 1984. Then you iterate through the cursor and grab the x & y from the shape, see here.

Answer 3

To elaborate on James's suggestion, here is a minimal code example using Python/arcpy:

import arcpy

def main():
    projectedPointFC = r'c:\point_test.shp'
    desc = arcpy.Describe(projectedPointFC)
    shapefieldname = desc.ShapeFieldName

    rows = arcpy.SearchCursor(projectedPointFC, r'', \
                              r'GEOGCS["GCS_WGS_1984",' + \
                              'DATUM["D_WGS_1984",' + \
                              'SPHEROID["WGS_1984",6378137,298.257223563]],' + \
                              'PRIMEM["Greenwich",0],' + \
                              'UNIT["Degree",0.017453292519943295]]')

    for row in rows:
        feat = row.getValue(shapefieldname)
        pnt = feat.getPart()
        print pnt.X, pnt.Y

if __name__ == '__main__':
    main()

Answer 4

Whether you call it projection or not, I am pretty sure that by definition, when you are translating the coordinate values from one spatial reference system to another, you are re/un-projecting.

I am not that familiar with ArcPy, but in arcgisscripting at 9.3, you would have to project the whole feature class.

Depending on how complex of a projection/transormation algorithm you need, you could always roll your own projection for the coordinates in basic python math. This would allow you to to coordinate value projection at the feature level.

If you were open to using the OGR python bindings, you can project at the feature level within something like a 'search cursor'.

Answer 5

At ArcPy 10.0 there is no ability to project individual geometries. However, you can create a feature set (or an in-memory feature class) and project that instead of a full-blown feature class in a workspace on disk or in a database somewhere.

Answer 6

The main reason I can see not wanting to create a feature class is because arcpy.CreateFeatureclass_management can be slow. You can also use arcpy.da.NumPyArrayTofeatureClass, which is more or less instant for in_memory feature classes:

In [1]: import arcpy

In [2]: import numpy as np

In [3]: geosr = arcpy.SpatialReference('Geographic Coordinate Systems/Spheroid-based/WGS 1984 Major Auxiliary Sphere')

In [4]: tosr = arcpy.SpatialReference('Projected Coordinate Systems/World/WGS 1984 Web Mercator (auxiliary sphere)')

In [5]: npai=list(enumerate(((-115.12799999956881, 36.11419999969922), (-117, 38.1141))))

In [6]: npai
Out[6]: [(0, (-115.12799999956881, 36.11419999969922)), (1, (-117, 38.1141))]

In [7]: npa=np.array(npai, np.dtype(([('idfield', np.int32), ('XY', np.float, 2)])))

In [8]: npa
Out[8]: 
array([(0, [-115.12799999956881, 36.11419999969922]),
       (1, [-117.0, 38.1141])], 
      dtype=[('idfield', '<i4'), ('XY', '<f8', (2,))])

In [9]: fcName = arcpy.CreateScratchName(workspace='in_memory', data_type='FeatureClass')

In [10]: arcpy.da.NumPyArrayToFeatureClass(npa, fcName, ['XY'], geosr)

In [11]: with arcpy.da.SearchCursor(fcName, 'SHAPE@XY', spatial_reference=tosr) as cur:
    ...:     print list(cur)
    ...:     
[((-12815990.336048, 4316346.515041453),), ((-13024380.422813002, 4595556.878958654),)]

Answer 7

import arcpy
dsc = arcpy.Describe(FC)
cursor = arcpy.UpdateCursor(FC, "", "Coordinate Systems\Geographic Coordinate   Systems\World\WGS 1984.prj")
for row in cursor:
  shape=row.getValue(dsc.shapeFieldName)
  geom = shape.getPart(0)
  x = geom.X
  y = geom.Y
  row.setValue('LONG_DD', x)
  row.setValue('LAT_DD', y)
  cursor.updateRow(row)

del cursor, row