Obtenir la latitude et la longitude du point projeté à l'aide d'ArcPy? [fermé]

13

J'ai une entité ponctuelle dans une classe d'entités à laquelle ArcPy accède. Le point est projeté mais je dois trouver un moyen efficace d'obtenir la latitude et la longitude non projetées pour ce point.

Existe-t-il une méthode autre que la reprojection (dé-projection), l'obtention d'un curseur de recherche sur la nouvelle classe d'entités, la recherche de l'entité, puis la suppression du lat / lon de la forme de l'entité?

Kenton W
la source

Réponses:

6

Le SearchCursor prend en charge la spécification d'une référence spatiale - dans ce cas, vous souhaitez un système de coordonnées géographiques, tel que WGS 1984. Ensuite, vous parcourez le curseur et saisissez les x et y de la forme, voir ici .

James
la source
6

La plupart des autres réponses ont toutes été publiées lorsque ArcGIS 10.0 était le dernier logiciel. Dans ArcGIS 10.1, de nombreuses nouvelles fonctionnalités ArcPy sont devenues disponibles. Cette réponse tire parti de cette nouvelle fonctionnalité. Il ne conviendra pas à 10.0 mais offre des performances et des fonctionnalités accrues pour 10.1 et versions ultérieures.

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

Cet extrait de code utilise le wkid pour créer un objet de référence spatiale plutôt que de taper une représentation sous forme de chaîne, utilise les curseurs d'accès aux données plus modernes et projette les objets de géométrie individuels à l'aide de la méthode projectAs () .

GeoSharp
la source
Bonne réponse. Je suggérerais simplement de changer X et Y dans l'impression, car dans WGS84 l'ordre commun est lat / long
radouxju
encore plus simple, faites-le. srs = arcpy.SpatialReference (4326) xy_coords = arcpy.da.FeatureClassToNumPyArray (input_feature_class, 'SHAPE @ XY', spatial_reference = srs) print (xy_coords)
dfresh22
5

Pour développer la suggestion de James, voici un exemple de code minimal utilisant 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
la source
4

Que vous l'appeliez projection ou non, je suis à peu près sûr que par définition, lorsque vous traduisez les valeurs de coordonnées d'un système de référence spatial à un autre, vous êtes en train de re / dé-projeter.

Je ne suis pas très familier avec ArcPy, mais dans arcgisscripting à 9.3, vous devrez projeter toute la classe d'entités.

Selon la complexité d'un algorithme de projection / transormation dont vous avez besoin, vous pouvez toujours rouler votre propre projection pour les coordonnées en mathématiques de base de python. Cela vous permettrait de coordonner la projection des valeurs au niveau de l'entité.

Si vous étiez ouvert à l'utilisation des liaisons Python OGR, vous pouvez projeter au niveau de la fonctionnalité dans quelque chose comme un «curseur de recherche».

DavidF
la source
Malheureusement, je ne peux pas utiliser des éléments non ESRI avec le script que j'utilise. Même si ESRI utilise OGR et GDAL (ne le dites à personne, n'est-ce pas?) ...
Kenton W
En fait, le meilleur itinéraire pourrait être de comprendre comment utiliser PROJ4 directement sur les coordonnées d'entrée d'une manière ou d'une autre.
Kenton W
@Kenton - Cela inclut-il également votre propre algorithme personnalisé (basé sur le code existant)? Si vous avez besoin de convertir UTM -> WGS84, j'ai du code pour le faire en python que je pourrais poster. Alternativement, vous pouvez extraire l'algorithme requis de Proj4 et l'utiliser à la place. Ou si vous êtes vraiment contraint d'utiliser du code ESRI (et que vous ne voulez pas projeter une classe d'entités entière comme suggéré), écrivez une bibliothèque C simple à projeter à l'aide d'ArcObjects, puis appelez-la depuis Python à l'aide de ctypes. Ou restez avec arcpy et projetez une classe d'entités entière :(
Sasa Ivetic
@Kenton - La recherche rapide renvoie pyproj ( code.google.com/p/pyproj ), vous pouvez y voir un exemple d'utilisation de python pour appeler la bibliothèque Proj4.
Sasa Ivetic
@Kenton - S'il s'agit d'une projection UTM NAD83 => WGS84 géographique sans transformation de datum, vous devriez être capable d'implémenter l'algorithme en python pur. Les équations sont dans le livre de Snyder: onlinepubs.er.usgs.gov/djvu/PP/PP_1395.pdf J'ai une fonction Oracle PL / SQL qui fait cela si vous voulez le code. J'avais l'intention de porter cette fonction sur Python, mais généralement, j'utilise simplement ogr / osr ...
DavidF
4

ArcPy 10.0 ne permet pas de projeter des géométries individuelles. Cependant, vous pouvez créer un ensemble d'entités (ou une classe d'entités en mémoire) et le projeter à la place d'une classe d'entités complète dans un espace de travail sur le disque ou dans une base de données quelque part.

Philippe
la source
c'est exactement ce que j'espérais éviter. Me fait souhaiter la puissance que vous pouvez obtenir dans .Net avec ArcObjects ...
Kenton W
0

La principale raison pour laquelle je vois que je ne souhaite pas créer une classe d'entités est que l'arcpy.CreateFeatureclass_management peut être lent. Vous pouvez également utiliser arcpy.da.NumPyArrayTofeatureClass, qui est plus ou moins instantané pour les classes d'entités 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
la source
-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
ThinkSpatially
la source