Extraindo coordenadas de vértices poligonais no ArcMap?

25

Eu tenho cerca de uma dúzia de polígonos em uma classe de recurso carregada no ArcMap 10, tudo no WGS 1984 geográfico.

Como obtenho facilmente as coordenadas associadas a cada vértice de cada polígono nessa classe de recurso?

Idealmente, eu gostaria que eles fossem tabulados em formato de planilha.

hpy
fonte

Respostas:

12

Isso funciona com uma licença padrão do ArcGIS:

desc = arcpy.Describe(fcl)
shapefieldname = desc.ShapeFieldName

gebieden = arcpy.UpdateCursor(fcl)

for gebied in gebieden:
    polygoon = gebied.getValue(shapefieldname)
    for punten in polygoon:
        for punt in punten:
            print punt.X, punt.Y
Bert H.
fonte
4

Esta é outra maneira de fazer isso usando um da.SearchCursor :

import arcpy
fc=r'C:\TEST.gdb\polygons123'

with arcpy.da.SearchCursor(fc,['OID@','SHAPE@']) as cursor:
    for row in cursor:
        array1=row[1].getPart()
        for vertice in range(row[1].pointCount):
            pnt=array1.getObject(0).getObject(vertice)
            print row[0],pnt.X,pnt.Y

Resultando em ObjectID, X e Y, que podem ser copiados para o Excel:

...
1 537505.894287 6731069.60889
1 537533.516296 6731078.20947
1 537555.316528 6731082.53589
1 537562.501892 6731085.47913
1 537589.395081 6731070.52991
1 537617.062683 6731058.29651
2 537379.569519 6729746.16272
2 537384.81311 6729746.06012
2 537396.085327 6729748.62311
2 537404.065674 6729752.75311
2 537425.145325 6729773.72931
2 537429.842102 6729777.07129
2 537442.971313 6729780.10651
2 537450.27533 6729780.51611
...
BERA
fonte
Pergunta rápida, o que está acontecendo na parte 'array1.getObject (0) .getObject (vértice)'?
Rex
@Rex, consulte a seção de ajuda do Array: pro.arcgis.com/en/pro-app/arcpy/classes/array.htm . A partir da matriz, busco
BERA
Ótimo, isso me ajuda a entender um pouco melhor. Então, por que existe uma matriz em uma matriz e não apenas uma matriz de pontos? Existe mais alguma coisa na primeira matriz? O que array1.getObject (1) forneceria a você?
Rex
Isso ocorre porque você obtém todas as partes em "row [1] .getPart ()"; portanto, sua primeira matriz são as diferentes partes de um recurso de várias partes. Então você só teria algo além de array1.getObject (0) se tivesse um recurso de várias partes?
Rex
3

Experimente as ferramentas de geo-assistente das tecnologias espaciais. Possui várias ferramentas gratuitas que podem fazer o que você deseja. Experimente as coordenadas do polígono get. Ou polígono para pontos

e assistentes geográficos

Brad Nesom
fonte
2

O script python a seguir (que requer o ArcGIS 10.1 ou posterior) utiliza arcpy.dapara pegar um shapefile como entrada e criar uma planilha com uma entrada para cada vértice em cada polígono presente no .shp (e acredito que funcione com licenças de arcgis de nível inferior) . Os IDs de objetos e sequências vinculam os pontos a uma posição específica em um polígono específico.

H / t para @PaulSmith nesta postagem: obtenha todos os pontos de uma polilinha para destacar a explode_to_pointsopção na arcpy.da.FeatureClassToNumPyArrayferramenta

import os
import csv
import arcpy
from os import path
from arcpy import da
from arcpy import env

env.overwriteOutput = True
env.workspace = '/folder/containing/your/shp/here'

polygon_shp = path.join(env.workspace, 'your_shapefile_name.shp')
vertex_csv_path = 'your/csv/path/here/poly_vertex.csv'

def getPolygonCoordinates(fc):
    """For each polygon geometry in a shapefile get the sequence number and
    and coordinates of each vertex and tie it to the OID of its corresponding
    polygon"""

    vtx_dict = {}
    s_fields = ['OID@', 'Shape@XY']
    pt_array = da.FeatureClassToNumPyArray(polygon_shp, s_fields, 
        explode_to_points=True)

    for oid, xy in pt_array:
        xy_tup = tuple(xy)
        if oid not in vtx_dict:
            vtx_dict[oid] = [xy_tup]
        # this clause ensures that the first/last point which is listed
        # twice only appears in the list once
        elif xy_tup not in vtx_dict[oid]:
            vtx_dict[oid].append(xy_tup)


    vtx_sheet = []
    for oid, vtx_list in vtx_dict.iteritems():
        for i, vtx in enumerate(vtx_list):
            vtx_sheet.append((oid, i, vtx[0], vtx[1]))

    writeVerticesToCsv(vtx_sheet)

def writeVerticesToCsv(vtx_sheet):
    """Write polygon vertex information to csv"""

    header = (
        'oid',          'sequence_id', 
        'x_coordinate', 'y_coordinate')

    with open(vertex_csv_path, 'wb') as vtx_csv:
        vtx_writer = csv.writer(vtx_csv)
        vtx_writer.writerow(header)

        for row in vtx_sheet:
            vtx_writer.writerow(row)

getPolygonCoordinates(polygon_shp)

Também escrevi um script que aborda especificamente os requisitos de: Inserir coordenadas de vértices no polígono, que é marcado como duplicado para esta pergunta, esse código está abaixo:

import os
import arcpy
from os import path
from arcpy import da
from arcpy import env
from arcpy import management

env.overwriteOutput = True
env.workspace = '/folder/containing/your/shp/here'

polygon_shp = path.join(env.workspace, 'your_shapefile_name.shp')
file_gdb = 'your/file/gdb/path/here/temp.gdb'

def addVerticesAsAttributes(fc):
    """Add the x,y coordinates of vertices as attributes to corresponding 
    features.  The coordinates will be in the order the appear in the geometry"""

    polygon_copy = createGdbFcCopy(fc)

    vtx_dict = {}
    s_fields = ['OID@', 'Shape@XY']
    pt_array = da.FeatureClassToNumPyArray(polygon_copy, s_fields, 
        explode_to_points=True)

    for oid, xy in pt_array:
        xy_tup = tuple(xy)
        if oid not in vtx_dict:
            vtx_dict[oid] = [xy_tup]
        # this clause ensures that the first/last point which is listed
        # twice only appears in the list once
        elif xy_tup not in vtx_dict[oid]:
            vtx_dict[oid].append(xy_tup)

    # find that largest number of points that exist within a polygon to determine 
    # the number of fields that need to be added to the shapefile
    max_vertices = 0
    for vtx_list in vtx_dict.values():
        if len(vtx_list) > max_vertices:
            max_vertices = len(vtx_list)

    xy_fields = addXyFields(polygon_copy, max_vertices)

    u_fields = ['OID@'] + xy_fields
    with da.UpdateCursor(polygon_copy, u_fields) as cursor:
        oid_ix = cursor.fields.index('OID@')
        for row in cursor:
            xy_ix = oid_ix + 1
            for vtx in vtx_dict[row[oid_ix]]:
                for coord in vtx:
                    row[xy_ix] = coord
                    xy_ix += 1

            cursor.updateRow(row)

def createGdbFcCopy(fc):
    """Create copy of the input shapefile as a file geodatabase feature class,
    because a huge number of fields may be added to the fc this preferable to shp"""

    if not arcpy.Exists(file_gdb):
        management.CreateFileGDB(path.dirname(file_gdb), 
            path.basename(file_gdb))

    polygon_copy = path.join(file_gdb, 'polygon_shp_copy')
    management.CopyFeatures(polygon_shp, polygon_copy)
    return polygon_copy

def addXyFields(fc, vtx_count):
    """Add fields to the feature class that will hold the x, y coordinates for each
    vertex, the number of fields is twice the number of most vertices in any polygon"""

    field_list = []
    f_type = 'DOUBLE'
    for i in range(1, vtx_count+1):
        f_names = ['x{0}'.format(i), 'y{0}'.format(i)]
        for fn in f_names:
            management.AddField(fc, fn, f_type)

        field_list.extend(f_names)

    return field_list

addVerticesAsAttributes(polygon_shp)
Grant Humphries
fonte
O primeiro script python faz o que o hpy pediu! Funciona muito rápido! Obrigado Grant Humphries!
ArisA
1

Ainda não concluí a solução, mas parece que você pode usar esta ferramenta:

Conversão> JSON> Recursos para JSON.

Isso converterá seu shapefile (no meu caso, 81 polígonos) em um arquivo JSON. Você pode abrir isso com um editor de texto para ver se todos os vértices estão listados para cada polígono.

Além disso, a biblioteca padrão do python (import json) trata os objetos json como dicionários. Você pode simplesmente percorrer seus dicionários para gravar os valores dos vértices (e quaisquer outros atributos desejados) em um arquivo csv. Se eu conseguir funcionar, voltarei e publicarei o soln.

nlb
fonte
0

Eu só precisava das coordenadas xey para a polilinha e o polígono. Eu usei o ToolBox -> Ferramentas de Gerenciamento de Dados -> Recursos -> Recurso a Ponto. Isso criou um arquivo de forma de ponto e, em seguida, usei as coordenadas XY adicionadas no mesmo menu Recursos para gerar coordenadas XY. Então eu extraí as informações da tabela de atributos shapes em uma planilha do excel. Isso resolveu o meu problema, não tenho certeza se você está procurando o mesmo.

Shashikiran
fonte
Isso não fornecerá apenas um ponto (X, Y) por polilinha / polígono quando a pergunta estiver pedindo um ponto (X, Y) para cada vértice em cada polilinha / polígono?
PolyGeo
-2

Aqui está uma solução alternativa em tempos desesperados:

  • Comece a editar a classe de recurso ou o shapefile
  • Selecione o recurso de polígono e clique com o botão direito do mouse em 'Editar vértices'
  • Clique com o botão direito do mouse em um dos vértices e selecione 'Propriedades do esboço'
  • Um menu desdobrável aparecerá com as coordenadas dos vértices listados
  • Tire uma captura de tela da lista de coordenadas
  • Cole a captura de tela no seu editor de fotos / fotos favorito e salve como jpeg / png / bmp etc.
  • Google 'OCR on-line gratuito' Escolha um dos resultados (alguns são melhores que outros)
  • Carregue seu arquivo da captura de tela coordenada e converta
  • Escolha o seu tipo de arquivo de saída (txt, Excel etc)
  • Verifique os resultados, pois alguns conversores de OCR são lixo !!!
  • Use os dados X, Y adicionados no Arcmap para criar um conjunto de dados de pontos.

Esse método é adequado para conjuntos de dados pequenos, mas a dependência / limitações dos conversores de OCR são a principal preocupação. Use com cuidado.

Mike
fonte