Calculando a orientação de cada lado do polígono usando o ArcPy?

9

Quero examinar a orientação de cada linha em um polígono para que eu possa calcular sua exposição solar. Cada polígono representa um edifício e tem uma altura associada. No momento, só quero considerar a orientação e, posteriormente, considerarei os problemas de sombreamento.

Uma abordagem que pensei foi dividir o polígono em linhas e calcular a orientação de cada linha, mas a dificuldade é que preciso identificar a face externa dessa linha. Embora a maioria dos polígonos sejam simples figuras quadridimensionais com linhas retas, há um pequeno número em que esse não é o caso (uma questão que eu apenas quero considerar, mas ainda não precisamos resolver).

Eu estou familiarizado com python e planejava fazer tudo isso a partir de um script.

arcpy arcgis-10.0 djq
fonte
11
As respostas a esta pergunta ajuda: gis.stackexchange.com/questions/1886/...
Sean
@ Sean - obrigado, sim, isso ajuda em geral. Não estou familiarizado com os comandos específicos do ArcGIS que podem ser usados ​​para fazer isso. Além disso, a questão de saber qual a face interna / externa da linha em um polígono ainda permanece.
DJQ
Quando você diz 'orientação', ainda não o definiu bem o suficiente para mim - qual é a orientação de um hexágono perfeito, por exemplo?
Dan S.
@ Dan S. Acabei de editar minha resposta para incluir a orientação de cada linha no polígono.
DJQ
desculpe respondedores - não consegui atribuir a recompensa corretamente ontem.
DJQ

Respostas:

6

Se você quer apenas orientação majoritária, confira a resposta de @Mapperz acima.

Caso contrário, como você disse, você poderá dividir os polígonos em linhas usando a ferramenta Polígono em linha . Isso adiciona um campo FID esquerdo e direito, onde o campo é -1 se não houver polígono externo - isso pode causar um pouco de confusão, se seus edifícios forem adjacentes ou se sobreporem.

A partir daí, você pode dividir suas linhas em todos os vértices (talvez use Dividir em linhas COGO ) e depois calcular os ângulos em cada uma das linhas (potencialmente atualizando os atributos COGO ).

Supondo que você tenha seu campo de ângulo calculado a partir do norte, o aspecto estará correto onde o left_FID é -1 e, para obter o aspecto quando o right_FID for -1, adicione 180 °. Depois, com base no FID original, você pode agregar, obter o aspecto majoritário com base no comprimento etc.

A ferramenta Polígono para linha é programável, (até onde eu sei), as ferramentas COGO não são, portanto você teria que criar algo lá.

Espero que isto ajude!

om_henners
fonte
2
Plug descarado e um comentário: primeiro, se você não possui uma licença Info, code.google.com/p/boundary-generator faz mais ou menos o mesmo que a ferramenta Polígono para Linha. Segundo - isso fornece muito mais informações do que apenas olhar para algum tipo de orientação geral do polígono geral ... por exemplo, você pode usar a tabela resultante para fazer cálculos muito mais precisos da exposição ao sol, levando em consideração o aspecto de cada parede individual .
Dan S.
5

Encontrando orientação

Em um script, o polígono estará disponível como um conjunto de anéis - um anel externo e zero ou mais anéis internos - com cada anel representado por uma tupla de vetores ciclicamente ordenada (v [0], v 1 , ... , v [m-1], v [m] = v [0]). Cada vetor fornece as coordenadas de um vértice (sem dois vértices sucessivos coincidentes). A partir disso, é simples, como outros já apontaram, obter vetores normais (isto é, vetores perpendiculares às direções das arestas):

n [i] = t (v [i + 1] - v [i]).

A operação "t" gira um vetor 90 graus no sentido anti - horário :

t ((x, y)) = (y, -x).

Somente as direções desses vetores normais são importantes, portanto, redimensione-as para ter comprimento unitário: um vetor (x, y) é redimensionado para (x / s, y / s) em que s = Sqrt (x ^ 2 + y ^ 2) (que é o comprimento da aresta correspondente). A partir de agora, vamos assumir que isso foi feito. Escreva os componentes dos vetores normais da unidade resultante como

n [i] = (u [i], v [i]), i = 0, 1, ..., m-1.

Discriminando do lado de fora

Como você observa, isso deixa uma ambiguidade direcional: devemos usar n [i] ou -n [i]? Qual aponta para o exterior? Esta pergunta é equivalente a encontrar o grau do mapa de Gauss . Para calcular, você precisa somar os ângulos pelos quais as direções normais mudam à medida que você marcha em torno de um anel. Como os vetores normais possuem comprimento unitário, o cosseno do ângulo entre duas arestas sucessivas é

Cos (q_i) = n [i]. n [i + 1] = u [i] * u [i + 1] + v [i] * v [i + 1], i = 0, 1, ..., m-1.

(Defina n [m] = n [0].)

O seno do ângulo entre duas arestas sucessivas é

Sin (q_i) = n [i]. t (n [i + 1]) = u [i] * v [i + 1] - v [i] * u [i + 1].

(Observe que esses cálculos exigem apenas somas, diferenças e produtos até o momento.) A aplicação da função tangente inversa principal (ATan2) a qualquer par (cosseno, seno) fornece o ângulo q_i entre -180 e 180 graus. A soma desses ângulos para i = 0, 1, ..., n-1 produz (até erro de ponto flutuante) a curvatura total do anel, que deve ser um múltiplo de 360 ​​graus; para um anel fechado sem interseção, ele será +360 ou -360. No primeiro caso, o grau é 1 e, no segundo caso, o grau é -1. As normais são todas orientadas para o exterior quando o grau do anel externo é +1 e os graus dos anéis internos são -1. Reordene-os anel a anel, conforme necessário, de acordo com esta regra. Ou seja, se o grau de qualquer anel for o oposto ao necessário, anule todas as normais desse anel. Agora você pode prosseguir com seus cálculos de insolação.

whuber
fonte
Obrigado por uma resposta muito abrangente. Quando você diz: 'Em um script, o polígono estará disponível como um conjunto de anéis'. Como o polígono está disponível? Embora familiarizado com alguns scripts python, não sei como interpretar o polígono dessa maneira. Estou lendo isso devagar e tentando entender, mas tenho dificuldade em traduzir algumas das explicações em pseudocódigo que posso escrever.
DJQ
Algumas notas sobre esta resposta: As APIs GIS típicas sempre apresentarão a parte externa de um polígono na ordem anti-horária, IIRC, o que significa que você não precisa fazer essa discriminação discriminatória de fora para dentro - use rotações no sentido horário nos anéis externos e anti-horário nos orifícios. Um esclarecimento para o celenius: o bit 'conjunto de anéis' é como a maioria das APIs, incluindo o python do ArcGIS, permite decompor um polígono em seus segmentos de linha constituintes - um anel é uma curva fechada deles. Desde polígonos pode apoiar ilhas e buracos, você pode ter vários anéis ...
Dan S.
@ Dan Gostaria que fosse o caso de GIS manter tais representações consistentes. Ao longo dos anos, o software ESRI flutuou entre os polígonos de orientação negativa e positiva, deixando-os sem orientação. Neste ponto, não tenho certeza se confiaria mesmo em declarações documentadas sobre sua abordagem atual: eu seria cauteloso. Custa pouco para verificar a orientação depois que você já processou todas as arestas em um anel.
whuber
.. bem, eu fiz esse pequeno aviso do IIRC, felizmente. ;) Não é tão frequente que eu faça mais coisas em nível de geometria na pilha ESRI.
Dan S.
@ Dan S. - Ainda estou um pouco incerto sobre como isso pode ser programado em python. Existem diretrizes para isso?
DJQ
1

Isso pode ajudar?

julien
fonte
Esse é um artigo interessante que você desenterrou. No entanto, nenhum dos métodos descritos lá é relevante para o cálculo da exposição solar (a menos que o cálculo pretenda ser uma aproximação bruta, o que não parece ser o caso aqui).
whuber
1

/ * Talvez isso ajude:

Azimute - pi / 2 é a orientação voltada para fora dos lados de um polígono de RHR:

Aqui está um exemplo do PostGIS, você pode criar a tabela bldg117862 usando a instrução no final. O SRID é EPSG 2271 (PA StatePlane North Feet) e a geometria é um multipolígono. Para visualizar no ArcGIS 10, cole consultas / subconsultas em uma conexão Query Layer para postgis após criar a tabela bldg117862. * /

- === INÍCIO DA CONSULTA ===

/ * A consulta externa fornece a orientação dos ortogonais externos e cria linhas ortogonais externas de comprimentos iguais aos dos lados do ponto médio dos lados.

As direções dominantes serão a soma do comprimento, agrupado por orientação, em ordem decrescente * /

SELECT line_id como side_id, comprimento, graus (ortoaz) como orientação, st_makeline (st_setsrid (st_line_interpolate_point (geom, .5), 2271), st_setsrid (st_makepoint (st_x (st_line_interpolate_point (geom, .5)) + (length * (sin (. orthoaz))), st_y (st_line_interpolate_point (geom, .5)) + (length * (cos (orthoaz)))), 2271)) como geom de

--a próxima subconsulta externa cria linhas a partir dos pares de pontos dos lados, calcula o azimute (ortoaz) da ortogonal externa para cada segmento

(SELECT bldg2009gid, line_id, st_length (st_makeline (ponto inicial, ponto final)) :: numeric (10,2) como comprimento, azimute (ponto inicial, ponto final), azimute (ponto inicial, ponto final) - pi () / 2 como ortoaz, st_makeline ( ponto inicial) e final) como geom de

/ * subconsulta mais interna - use generate_series () para decompor polígonos de construção em pares de pontos de extremidade / ponto final dos lados - note1 - force a regra da mão direita para garantir a orientação comum de todos os lados do polígono note2 - exemplo usa multipolígono, para polígono the geometryn () pode ser removido */

(SELECT generate_series (1, npoints (exteriorring (geometryn (st_forceRHR (geom), 1))) - 1) como line_id, gid como bldg2009gid, pointn (exteriorring (geometryn (st_forceRHR (geom), 1)), generate_series (1, npoints (exteriorring (geometryn (st_forceRHR (geom), 1))) - 1)) como ponto inicial, pointn (exteriorring (geometryn (st_forceRHR (geom), 1)), generate_series (2, npoints (exteriorring (geometryn (st_forceRHR (geom), 1)) ), 1))))) como ponto final de bldg117862) como t1) como t2

- === FIM DA CONSULTA ===

- as instruções de criação / inserção da tabela bldg117862

Defina STANDARD_CONFORMING_STRINGS para ON; SELECT DropGeometryColumn ('', 'bldg117862', 'geom'); TABELA DE QUEDA "bldg117862"; INÍCIO; CREATE TABLE "bldg117862" (gid serial PRIMARY KEY, "motherpin" varchar (14), "taxpin" varchar (14), "taxpin" varchar (14), "status" varchar (15), "area" numeric, "prev_area" numeric, "pct_change" numeric, "imagem" varchar (133), "página de mapa" varchar (6), "sref_gid" int4, "endereço_e" varchar (19), "a_address" varchar (19), "perim" numérico, "cartão" int4, "a_addnum" int4, "e_street" varchar (50), "a_street" varchar (50), "e_hsnum" varchar (10)); SELECT AddGeometryColumn ('', 'bldg117862', 'geom', '2271', 'MULTIPOLYGON', 2); 0106000020DF080000010000000103000020DF080000010000000B0000008C721D6C98AC34415E2C5BB9D3E32541AE56DE17BEAC34410613E5A0A0E325411AB6C794AEAC3441BA392FE372E32541C89C38429DAC3441643857628AE325418C299A9095AC3441F66C29B573E32541983F02087EAC34413080AA9F93E325419BAC3C0A86AC3441AC1F3B3DABE32541803A40B974AC3441E8CF3DB9C2E325413E3758C186AC3441D0AAB0E7F7E325410AAAA5429BAC3441BA971217DCE325418C721D6C98AC34415E2C5BB9D3E32541' ); CREATE INDEX "bldg117862_geom_gist" ON "bldg117862" usando gist ("geom" gist_geometry_ops); FIM;


fonte
0

Supondo que a orientação dos segmentos de linha seja constante em um polígono, você pode calcular o rumo (título) de um vetor perpendicular a cada segmento de linha. Agora não tenho tempo para me basear no código, mas se você precisar de matemática, ela pode ser facilmente fornecida :-)

WolfOdrade
fonte