Получение широты и долготы проектируемой точки с помощью ArcPy? [закрыто]

13

У меня есть точечный объект в классе объектов, к которому обращается ArcPy. Точка проецируется, но мне нужно найти эффективные средства для получения непрогнозируемой широты и долготы для этой точки.

Существует ли какой-либо метод, кроме перепроецирования (отмены проецирования), получения курсора поиска для нового класса пространственных объектов, поиска объекта, а затем получения широты / долготы от формы объекта?

Кентон W
источник

Ответы:

6

SearchCursor поддерживает указание пространственной привязки - в этом случае вам понадобится географическая система координат, такая как WGS 1984. Затем вы перебираете курсор и получаете x & y из фигуры, см. Здесь .

Джеймс
источник
6

Большинство других ответов были опубликованы, когда ArcGIS 10.0 была последней версией программного обеспечения. В ArcGIS 10.1 появилось много новых функций ArcPy. Этот ответ использует преимущества этой новой функциональности. Он не подходит для 10.0, но предлагает повышенную производительность и функциональность для 10.1 и более поздних версий.

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

Этот фрагмент кода использует wkid для создания объекта пространственной привязки, а не для ввода строкового представления, использует более современные курсоры доступа к данным и проецирует отдельные геометрические объекты с помощью метода projectAs () .

GeoSharp
источник
хороший ответ. Я бы просто предложил поменять X и Y на отпечатке, потому что в WGS84 общий порядок
широта
еще проще, просто сделай это. srs = arcpy.SpatialReference (4326) xy_coords = arcpy.da.FeatureClassToNumPyArray (input_feature_class, 'SHAPE @ XY' ,atial_reference = srs) print (xy_coords)
dfresh22
5

Чтобы развить предложение Джеймса, вот минимальный пример кода с использованием 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()
Аллан Адэйр
источник
4

Называете ли вы это проекцией или нет, я вполне уверен, что по определению, когда вы переводите значения координат из одной системы пространственной привязки в другую, вы повторно / не проектируете.

Я не очень знаком с ArcPy, но в arcgisscripting в 9.3 вам придется спроектировать весь класс объектов.

В зависимости от того, насколько сложен алгоритм проекции / преобразования, вы всегда можете свернуть свою собственную проекцию для координат в базовой математике Python. Это позволит вам координировать проекцию значения на уровне объекта.

Если вы были открыты для использования привязок Python OGR, вы можете проецировать на уровне объектов в нечто вроде «поискового курсора».

DavidF
источник
К сожалению, я не могу использовать не-ESRI вещи со сценарием, который я использую. Даже если даже ESRI использует OGR и GDAL (не говорите никому, верно?) ...
Кентон W
На самом деле, лучшим способом было бы выяснить, как использовать PROJ4 непосредственно для входных координат.
Кентон W
@Kenton - это также включает в себя ваш собственный алгоритм (на основе существующего кода)? Если вам нужно конвертировать UTM -> WGS84, у меня есть код для этого на python, который я мог бы опубликовать. Кроме того, вы можете извлечь требуемый алгоритм из Proj4 и использовать его вместо этого. Или, если вы действительно ограничены использованием кода ESRI (и не хотите проектировать весь класс объектов, как предложено), напишите простую библиотеку C для проектирования с использованием ArcObjects, а затем вызовите ее из Python с использованием ctypes. Или придерживайтесь arcpy и спроектируйте весь класс объектов :(
Sasa Ivetic
@Kenton - Быстрый поиск возвращает pyproj ( code.google.com/p/pyproj ), вы можете посмотреть пример использования python для вызова библиотеки Proj4.
Саша Иветик
@Kenton - если это UTM NAD83 => географическая проекция WGS84 без преобразования данных, вы должны иметь возможность реализовать алгоритм на чистом питоне. Уравнения приведены в книге Снайдера: onlinepubs.er.usgs.gov/djvu/PP/PP_1395.pdf У меня есть функция Oracle PL / SQL, которая делает это, если вам нужен код. Я хотел перенести эту функцию на Python, но обычно просто использую ogr / osr ...
DavidF
4

В ArcPy 10.0 нет возможности проецировать отдельные геометрии. Однако вы можете создать набор объектов (или класс объектов в памяти) и спроецировать его вместо полноценного класса объектов в рабочем пространстве на диске или где-нибудь в базе данных.

Филипп
источник
это именно то, чего я надеялся избежать. Заставляет меня пожелать силы, которую вы можете получить в .Net с ArcObjects ...
Кентон W
0

Основная причина, по которой я не хочу создавать класс пространственных объектов, заключается в том, что arcpy.CreateFeatureclass_management может работать медленно. Вы также можете использовать arcpy.da.NumPyArrayTofeatureClass, который более или менее мгновенен для классов объектов 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
источник
-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
источник