Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question

6 Answers 6

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. http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/SearchCursor/000v00000039000000/

share|improve this answer
    
+1, this is the way to do it using arcpy, specifically. –  Allan Adair Sep 21 '11 at 13:54

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'.

share|improve this answer
    
Unfortunately I can't use non-ESRI stuff with the script I'm using. Even though even ESRI uses OGR and GDAL (don't tell anyone, right?)... –  Kenton W Apr 6 '11 at 5:18
    
Actually, the better route might be to figure out how to use PROJ4 directly on the input coordinates somehow. –  Kenton W Apr 6 '11 at 5:19
    
@Kenton - Does that also include your own custom algorithm (based on existing code)? If you need to convert UTM -> WGS84, I have code to do that in python I could post. Alternately, you could extract the required algorithm from Proj4 and use that instead. Or if you're really constrained to using ESRI code (and don't want to project an entire feature class like suggested), write a simple C library to project using ArcObjects, then call it from Python using ctypes. Or stick with arcpy and project an entire feature class :( –  Sasa Ivetic Apr 6 '11 at 13:29
    
@Kenton - Quick search returns pyproj (code.google.com/p/pyproj), you could look at that for an example of how to use python to call the Proj4 library. –  Sasa Ivetic Apr 6 '11 at 13:32
    
@Kenton - If it is a UTM NAD83 => geographic WGS84 projection with no datum transform, you should be able to implement the algorithm in pure python. The equations are in Snyder's book: onlinepubs.er.usgs.gov/djvu/PP/PP_1395.pdf I have an Oracle PL/SQL function that does this if you want the code. I have been meaning to port this function to Python, but usually just use ogr/osr... –  DavidF Apr 6 '11 at 16:01

At ArcPy 10.0 there's 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.

share|improve this answer
    
which is exactly what I was hoping to avoid. Makes me wish for the power you can get in .Net with ArcObjects... –  Kenton W Apr 6 '11 at 5:17

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()
share|improve this answer

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),)]
share|improve this answer
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
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.