Uzyskujesz szerokość i długość geograficzną rzutowanego punktu za pomocą ArcPy? [Zamknięte]

13

Mam funkcję punktową w klasie obiektów, do której dostęp ma ArcPy. Punkt jest przewidziany, ale muszę znaleźć skuteczny sposób na uzyskanie nie rzuconej szerokości i długości geograficznej dla tego punktu.

Czy istnieje metoda inna niż ponowne rzutowanie (odznaczanie), uzyskiwanie kursora wyszukiwania w nowej klasie elementów, znajdowanie elementu, a następnie usuwanie kształtu lat / lon z kształtu elementu?

Kenton W.
źródło

Odpowiedzi:

6

SearchCursor obsługuje określenie odniesienia przestrzennego - w tym przypadku potrzebujesz Geograficznego Układu Współrzędnych, takiego jak WGS 1984. Następnie iterujesz kursor i pobierasz x i y z kształtu, patrz tutaj .

James
źródło
6

Większość pozostałych odpowiedzi została opublikowana, gdy ArcGIS 10.0 był najnowszym oprogramowaniem. W ArcGIS 10.1 udostępniono wiele nowych funkcji ArcPy. Ta odpowiedź wykorzystuje tę nową funkcjonalność. Nie będzie odpowiedni dla wersji 10.0, ale oferuje zwiększoną wydajność i funkcjonalność dla wersji 10.1 i nowszych.

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

Ten fragment kodu używa wkid do utworzenia przestrzennego obiektu odniesienia zamiast wpisywania ciągów znaków, używa bardziej nowoczesnych kursorów dostępu do danych i wyświetla poszczególne obiekty geometrii za pomocą metody projectAs () .

GeoSharp
źródło
niezła odpowiedź. Proponuję po prostu zamienić X i Y na wydruku, ponieważ w WGS84 wspólnym porządkiem jest lat / long
radouxju
jeszcze prościej, po prostu zrób to. srs = arcpy.SpatialReference (4326) xy_coords = arcpy.da.FeatureClassToNumPyArray (klasa_wejściowa, „SHAPE @ XY”, odniesienie_przestrzenne = srs) print (xy_coords)
dfresh22
5

Aby rozwinąć sugestię Jamesa, oto minimalny przykład kodu wykorzystujący 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
źródło
4

Bez względu na to, czy nazywasz to rzutowaniem, czy nie, jestem całkiem pewien, że z definicji, kiedy tłumaczysz wartości współrzędnych z jednego układu odniesienia przestrzennego na inny, ponownie / odrzucisz.

ArcPy nie znam się tak dobrze, ale w arcgisscripting w wersji 9.3 musiałbyś rzutować całą klasę funkcji.

W zależności od złożoności algorytmu projekcji / transormacji potrzebujesz zawsze rzucić własną projekcję dla współrzędnych w podstawowej matematyce python. Umożliwiłoby to koordynację projekcji wartości na poziomie obiektu.

Jeśli byłeś otwarty na używanie powiązań pytona OGR, możesz wyświetlać na poziomie funkcji w coś w rodzaju „kursora wyszukiwania”.

DavidF
źródło
Niestety nie mogę używać rzeczy innych niż ESRI ze skryptem, którego używam. Nawet jeśli nawet ESRI używa OGR i GDAL (nie mów nikomu, prawda?) ...
Kenton W
W rzeczywistości lepszym sposobem może być wymyślenie, w jaki sposób używać PROJ4 bezpośrednio na współrzędnych wejściowych.
Kenton W
@Kenton - Czy obejmuje to również Twój własny algorytm (oparty na istniejącym kodzie)? Jeśli potrzebujesz przekonwertować UTM -> WGS84, mam kod, aby to zrobić w pythonie, który mógłbym opublikować. Alternatywnie możesz wyodrębnić wymagany algorytm z Proj4 i użyć go zamiast tego. Lub jeśli jesteś naprawdę ograniczony do używania kodu ESRI (i nie chcesz rzutować całej klasy obiektów, jak sugerowane), napisz prostą bibliotekę C do projektu za pomocą ArcObjects, a następnie wywołaj ją z Pythona przy użyciu ctypów. Lub trzymaj się arcpy i wyświetlaj całą klasę funkcji :(
Sasa Ivetic
@Kenton - Szybkie wyszukiwanie zwraca pyproj ( code.google.com/p/pyproj ), możesz na to zobaczyć przykład użycia Pythona do wywoływania biblioteki Proj4.
Sasa Ivetic
@Kenton - Jeśli jest to projekcja UTM NAD83 => geograficzna WGS84 bez transformacji odniesienia, powinieneś być w stanie zaimplementować algorytm w czystym pythonie. Równania są w książce Snydera: onlinepubs.er.usgs.gov/djvu/PP/PP_1395.pdf Mam funkcję Oracle PL / SQL, która robi to, jeśli chcesz kodu. Chciałem przenieść tę funkcję do Pythona, ale zwykle po prostu używam ogr / osr ...
DavidF
4

W ArcPy 10.0 nie ma możliwości wyświetlania indywidualnych geometrii. Możesz jednak utworzyć zestaw funkcji (lub klasę funkcji w pamięci) i wyświetlić ją zamiast w pełni funkcjonalnej klasy funkcji w obszarze roboczym na dysku lub w bazie danych.

Philip
źródło
dokładnie tego chciałem uniknąć. Sprawia, że ​​życzę sobie mocy, którą można uzyskać. Net dzięki ArcObjects ...
Kenton W
0

Głównym powodem, dla którego nie widzę potrzeby tworzenia klasy elementów, jest to, że arcpy.CreateFeatureclass_management może być powolny. Możesz także użyć arcpy.da.NumPyArrayTofeatureClass, która jest mniej więcej natychmiastowa dla klas obiektów 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
źródło
-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
Pomyśl przestrzennie
źródło