¿Obtener la latitud y la longitud del punto proyectado usando ArcPy? [cerrado]

13

Tengo una entidad de puntos en una clase de entidad a la que ArcPy está accediendo. El punto está proyectado pero necesito encontrar un medio eficiente para obtener la latitud y longitud no proyectadas para este punto.

¿Existe otro método que no sea la reproyección (desproyección), obtener un cursor de búsqueda en la nueva clase de entidad, encontrar la entidad y luego obtener el lat / lon fuera de la forma de la entidad?

Kenton W
fuente

Respuestas:

6

SearchCursor admite la especificación de una referencia espacial; en este caso, desearía un Sistema de coordenadas geográficas, como WGS 1984. Luego itera a través del cursor y toma el x & y de la forma, vea aquí .

James
fuente
6

La mayoría de las otras respuestas se publicaron cuando ArcGIS 10.0 era el último software. En ArcGIS 10.1, mucha nueva funcionalidad de ArcPy estuvo disponible. Esta respuesta aprovecha esa nueva funcionalidad. No será adecuado para 10.0, pero ofrece un mayor rendimiento y funcionalidad para 10.1 y posteriores.

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

Este fragmento de código usa el wkid para crear un objeto de referencia espacial en lugar de escribir una representación de cadena, usa los cursores de acceso a datos más modernos y proyecta los objetos de geometría individuales usando el método projectAs () .

GeoSharp
fuente
buena respuesta. Simplemente sugeriría cambiar X e Y en la impresión, porque en WGS84 el orden común es lat / long
radouxju
aún más simple, solo haz esto. srs = arcpy.SpatialReference (4326) xy_coords = arcpy.da.FeatureClassToNumPyArray (input_feature_class, 'SHAPE @ XY', spatial_reference = srs) print (xy_coords)
dfresh22
5

Para ampliar la sugerencia de James, aquí hay un ejemplo de código mínimo usando 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()
Allan Adair
fuente
4

Ya sea que lo llames proyección o no, estoy bastante seguro de que, por definición, cuando estás traduciendo los valores de coordenadas de un sistema de referencia espacial a otro, estás re-proyectando.

No estoy tan familiarizado con ArcPy, pero en arcgisscripting en 9.3, tendría que proyectar toda la clase de entidad.

Dependiendo de qué tan complejo sea el algoritmo de proyección / transformación que necesite, siempre puede obtener su propia proyección para las coordenadas en matemáticas básicas de Python. Esto le permitiría coordinar la proyección del valor a nivel de entidad.

Si estaba abierto a usar los enlaces de Python de OGR, puede proyectar en el nivel de característica dentro de algo como un 'cursor de búsqueda'.

DavidF
fuente
Desafortunadamente no puedo usar cosas que no sean ESRI con el script que estoy usando. Aunque incluso ESRI usa OGR y GDAL (no se lo digas a nadie, ¿verdad?) ...
Kenton W
En realidad, la mejor ruta podría ser descubrir cómo usar PROJ4 directamente en las coordenadas de entrada de alguna manera.
Kenton W
@Kenton: ¿eso también incluye su propio algoritmo personalizado (basado en el código existente)? Si necesita convertir UTM -> WGS84, tengo un código para hacerlo en Python que podría publicar. Alternativamente, puede extraer el algoritmo requerido de Proj4 y usarlo en su lugar. O si realmente está limitado a usar el código ESRI (y no desea proyectar una clase de entidad completa como se sugiere), escriba una biblioteca C simple para proyectar usando ArcObjects, luego llámela desde Python usando ctypes. O quédese con arcpy y proyecte una clase de entidad completa :(
Sasa Ivetic
@Kenton: la búsqueda rápida devuelve pyproj ( code.google.com/p/pyproj ), puede ver eso como un ejemplo de cómo usar python para llamar a la biblioteca Proj4.
Sasa Ivetic
@Kenton: si se trata de una proyección UTM NAD83 => geográfica WGS84 sin transformación de datos, debería poder implementar el algoritmo en Python puro. Las ecuaciones están en el libro de Snyder: onlinepubs.er.usgs.gov/djvu/PP/PP_1395.pdf Tengo una función Oracle PL / SQL que hace esto si desea el código. He tenido la intención de portar esta función a Python, pero generalmente solo uso ogr / osr ...
DavidF
4

En ArcPy 10.0 no hay capacidad de proyectar geometrías individuales. Sin embargo, puede crear un conjunto de características (o una clase de entidad en memoria) y proyectarlo en lugar de una clase de entidad completa en un espacio de trabajo en el disco o en una base de datos en algún lugar.

Philip
fuente
que es exactamente lo que esperaba evitar. Me hace desear el poder que puede obtener .Net con ArcObjects ...
Kenton W
0

La razón principal por la que puedo ver que no quiero crear una clase de entidad es porque arcpy.CreateFeatureclass_management puede ser lento. También puede usar arcpy.da.NumPyArrayTofeatureClass, que es más o menos instantáneo para las clases de entidad in_memory:

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),)]
cwa
fuente
-1
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
Piense espacialmente
fuente