Abrufen von Breiten- und Längengrad des projizierten Punkts mit ArcPy? [geschlossen]

13

Ich habe ein Punkt-Feature in einer Feature-Class, auf die ArcPy zugreift. Der Punkt ist projiziert, aber ich muss ein effizientes Mittel finden, um den nicht projizierten Breiten- und Längengrad für diesen Punkt zu ermitteln.

Gibt es eine andere Methode als die Neuprojektion (Projektion aufheben), einen Suchcursor für die neue Feature-Class abzurufen, das Feature zu finden und das Lat / Lon von der Form des Features zu entfernen?

Kenton W
quelle

Antworten:

6

Der SearchCursor unterstützt die Angabe eines Raumbezugs. In diesem Fall möchten Sie ein geografisches Koordinatensystem, z. B. WGS 1984, verwenden. Durchlaufen Sie dann den Cursor und nehmen Sie das x & y aus der Form, siehe hier .

James
quelle
6

Die meisten anderen Antworten wurden veröffentlicht, als ArcGIS 10.0 die neueste Software war. Ab ArcGIS 10.1 sind viele neue ArcPy-Funktionen verfügbar. Diese Antwort nutzt diese neue Funktionalität. Es ist nicht für 10.0 geeignet, bietet jedoch eine höhere Leistung und Funktionalität für 10.1 und höher.

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

Dieser Code - Schnipsel , die verwendet wkid um ein räumliches Referenzobjekt zu erstellen , anstatt die Eingabe aus einer String - Darstellung, verwendet die modernere Datenzugriff Cursor und Projekte der einzelnen Geometrieobjekte mit Hilfe der projectAs () -Methode .

GeoSharp
quelle
gute Antwort. Ich würde einfach vorschlagen, X und Y im Druck zu
vertauschen
Noch einfacher, machen Sie es einfach. SRS = arcpy.SpatialReference (4326) xy_coords = arcpy.da.FeatureClassToNumPyArray (input_feature_class, '-Form @ XY', spatial_reference = SRS) print (xy_coords)
dfresh22
5

Um auf James 'Vorschlag näher einzugehen, hier ein minimales Codebeispiel mit 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
quelle
4

Ganz gleich, ob Sie es Projektion nennen oder nicht, ich bin mir ziemlich sicher, dass Sie per Definition beim Übersetzen der Koordinatenwerte von einem Raumbezugssystem in ein anderes eine Neu- / Neuprojektion durchführen.

Ich bin nicht so vertraut mit ArcPy, aber in ArcGisscripting in Version 9.3 müssten Sie die gesamte Feature-Class projizieren.

Je nachdem, wie komplex ein Projektions- / Transformationsalgorithmus ist, können Sie in der einfachen Python-Mathematik jederzeit eine eigene Projektion für die Koordinaten ausführen. Auf diese Weise können Sie die Wertprojektion auf Feature-Ebene koordinieren.

Wenn Sie offen für die Verwendung der OGR-Python-Bindungen waren, können Sie auf Feature-Ebene in so etwas wie einem "Suchcursor" projizieren.

DavidF
quelle
Leider kann ich mit dem Skript, das ich verwende, keine Nicht-ESRI-Inhalte verwenden. Obwohl auch ESRI OGR und GDAL verwendet (sagen Sie es niemandem, oder?) ...
Kenton W
Tatsächlich könnte die bessere Route darin bestehen, herauszufinden, wie PROJ4 direkt auf den Eingabekoordinaten verwendet wird.
Kenton W
@Kenton - Enthält das auch Ihren eigenen benutzerdefinierten Algorithmus (basierend auf vorhandenem Code)? Wenn Sie UTM konvertieren müssen -> WGS84, habe ich Code, um das in Python zu tun, das ich posten konnte. Alternativ können Sie den erforderlichen Algorithmus aus Proj4 extrahieren und diesen stattdessen verwenden. Oder wenn Sie sich wirklich darauf beschränken, ESRI-Code zu verwenden (und nicht die gesamte Feature-Class wie vorgeschlagen projizieren möchten), schreiben Sie eine einfache C-Bibliothek, um sie mit ArcObjects zu projizieren, und rufen Sie sie dann in Python mit ctypes auf. Oder bleiben Sie bei arcpy und projizieren Sie eine ganze Feature-Class :(
Sasa Ivetic
@Kenton - Die Schnellsuche gibt pyproj zurück ( code.google.com/p/pyproj ). Hier finden Sie ein Beispiel für den Aufruf der Proj4-Bibliothek mit Python.
Sasa Ivetic
@ Kenton - Wenn es sich um eine UTM NAD83 => geografische WGS84-Projektion ohne Datumstransformation handelt, sollten Sie in der Lage sein, den Algorithmus in reinem Python zu implementieren. Die Gleichungen befinden sich in Snyders Buch: onlinepubs.er.usgs.gov/djvu/PP/PP_1395.pdf Ich habe eine Oracle PL / SQL-Funktion, die dies ausführt , wenn Sie den Code möchten. Ich wollte diese Funktion nach Python portieren, benutze aber normalerweise nur ogr / osr ...
DavidF
4

In ArcPy 10.0 können keine einzelnen Geometrien projiziert werden. Sie können jedoch ein Feature-Set (oder eine speicherinterne Feature-Class) erstellen und dieses anstelle einer vollständigen Feature-Class in einem Arbeitsbereich auf der Festplatte oder in einer Datenbank irgendwo projizieren.

Philip
quelle
Genau das wollte ich vermeiden. Ich wünsche mir die Leistung, die Sie in .Net mit ArcObjects erhalten können ...
Kenton W
0

Der Hauptgrund, warum ich sehe, dass ich keine Feature-Class erstellen möchte, ist, dass arcpy.CreateFeatureclass_management langsam sein kann. Sie können auch arcpy.da.NumPyArrayTofeatureClass verwenden, was für in_memory-Feature-Classes mehr oder weniger augenblicklich ist:

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