Cálculo de la orientación de cada lado del polígono usando ArcPy?

9

Quiero examinar la orientación de cada línea en un polígono para poder calcular su exposición solar. Cada polígono representa un edificio y tiene una altura asociada. Por el momento, solo quiero considerar la orientación y luego consideraré los problemas de sombreado.

Un enfoque que pensé fue dividir el polígono en líneas y calcular la orientación de cada línea, pero la dificultad es que luego tengo que identificar la cara externa de esa línea. Aunque la mayoría de los polígonos son figuras simples de cuatro lados con líneas rectas, hay un pequeño número donde este no es el caso (un problema que solo quiero considerar, pero que aún no tengo que resolver).

Estoy familiarizado con Python y planeé hacer todo esto desde un script.

arcpy arcgis-10.0 djq
fuente
1
Las respuestas a esta pregunta ayudan: gis.stackexchange.com/questions/1886/…
Sean
@Sean - gracias, sí, eso ayuda en general. No estoy familiarizado con los comandos específicos de ArcGIS que se pueden usar para hacer esto. Además, el problema de saber cuál es la cara interior / exterior de la línea en un polígono aún permanece.
DJ
Cuando dices 'orientación', no lo has definido lo suficientemente bien para mí --- ¿cuál es la orientación de un hexágono perfecto, por ejemplo?
Dan S.
@Dan S. Acabo de editar mi respuesta para incluir la orientación de cada línea en el polígono.
djq
lo siento, respondedores: ayer no pude asignar la recompensa correctamente.
djq

Respuestas:

6

Si solo desea una orientación mayoritaria, consulte la respuesta de @Mapperz arriba.

De lo contrario, como usted dice, podría dividir los polígonos en líneas con la herramienta Polígono a línea . Esto agrega un campo FID izquierdo y derecho, donde el campo es -1 si no hay un polígono externo; sin embargo, esto puede causar un poco de suciedad si sus edificios son adyacentes o se superponen.

A partir de ahí, puede dividir sus líneas en cada vértice (tal vez use Dividir en líneas COGO ) y luego calcular los ángulos en cada una de las líneas (potencialmente actualizando los atributos COGO ).

Suponiendo que tiene su campo de ángulo calculado desde el Norte, el aspecto será correcto donde left_FID es -1, y para obtener el aspecto cuando right_FID es -1 solo agregue 180 °. Luego, en función del FID original, puede agregar, obtener el aspecto mayoritario en función de la longitud, etc.

La herramienta Polígono a línea es programable, (hasta donde yo sé) las herramientas COGO no lo son, por lo que tendría que encontrar algo allí usted mismo.

¡Espero que esto ayude!

om_henners
fuente
2
Plug desvergonzado, y un comentario: Primero, si no tiene una licencia de información, code.google.com/p/boundary-generator hace más o menos lo mismo que la herramienta Polígono a línea. Segundo: esto le brinda mucha más información que solo mirar algún tipo de orientación general del polígono general ... por ejemplo, puede usar la tabla resultante para hacer cálculos de exposición al sol mucho más precisos, teniendo en cuenta el aspecto de cada pared individual .
Dan S.
5

Encontrar orientación

Dentro de un script, el polígono estará disponible como un conjunto de anillos (un anillo externo y cero o más anillos internos) con cada anillo representado por una tupla de vectores ordenada cíclicamente (v [0], v 1 , ... , v [m-1], v [m] = v [0]). Cada vector da las coordenadas de un vértice (sin que coincidan dos vértices sucesivos). A partir de esto, es sencillo, como han señalado otros, obtener vectores normales (es decir, vectores perpendiculares a las direcciones de los bordes):

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

La operación "t" gira un vector 90 grados en sentido antihorario :

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

Solo importan las direcciones de estos vectores normales, así que reescalarlos para que tengan una longitud unitaria: un vector (x, y) reescala a (x / s, y / s) donde s = Sqrt (x ^ 2 + y ^ 2) (que es la longitud de su borde correspondiente). De ahora en adelante, supongamos que esto se ha hecho. Escriba los componentes de los vectores normales de la unidad resultante como

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

Discriminar desde afuera hacia adentro

Como comentan, esto deja una ambigüedad direccional: ¿deberíamos usar n [i] o -n [i]? ¿Cuál apunta hacia afuera? Esta pregunta es equivalente a encontrar el grado del mapa de Gauss . Para calcularlo, debes sumar los ángulos por los que cambian las direcciones normales a medida que marchas alrededor de un anillo. Debido a que los vectores normales tienen una unidad de longitud, el coseno del ángulo entre dos bordes sucesivos es

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].)

El seno del ángulo entre dos bordes sucesivos es

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

(Tenga en cuenta que estos cálculos solo requieren sumas, diferencias y productos hasta ahora). La aplicación de la función de tangente inversa principal (ATan2) a cualquier par (coseno, seno) da el ángulo q_i entre -180 y 180 grados. Sumando estos ángulos para i = 0, 1, ..., n-1 produce (hasta error de coma flotante) la curvatura total del anillo, que debe ser un múltiplo de 360 ​​grados; para un anillo cerrado sin intersección, será +360 o -360. En el primer caso el grado es 1 y en el segundo caso el grado es -1. Las normales están todas orientadas hacia afuera cuando el grado del anillo externo es +1 y los grados de los anillos internos son -1. Reoriéntelos anillo por anillo, según sea necesario, de acuerdo con esta regla. Es decir, si el grado de cualquier anillo es el opuesto al que se necesita, niegue todas las normales para ese anillo. Ahora puede continuar con sus cálculos de insolación.

whuber
fuente
Gracias por una respuesta muy completa. Cuando dices, 'Dentro de un script, el polígono estará disponible como un conjunto de anillos' ¿cómo está disponible el polígono? Aunque estoy familiarizado con algunas secuencias de comandos de Python, no sé cómo interpretar el polígono de esta manera. Estoy leyendo esto lentamente y tratando de entenderlo, pero tengo dificultades para traducir parte de la explicación en pseudocódigo que puedo escribir.
djq
Algunas notas sobre esta respuesta: las API de SIG típicas siempre presentarán el exterior de un polígono en el sentido contrario a las agujas del reloj, IIRC, lo que significa que no tiene que hacer esa elegante distinción desde afuera hacia adentro, solo use rotaciones en sentido horario en los anillos exteriores y en sentido antihorario en los agujeros. Una aclaración para celenius: el bit del 'conjunto de anillos' es cómo la mayoría de las API, incluida la pitón de ArcGIS, le permiten descomponer un polígono en sus segmentos de línea constituyentes: un anillo es una curva cerrada de ellos. Dado que los polígonos pueden soportar islas y agujeros, podría tener múltiples anillos ...
Dan S.
@Dan, desearía que los SIG mantuvieran representaciones tan consistentes. A lo largo de los años, el software ESRI ha fluctuado de un lado a otro entre polígonos orientados negativa y positivamente y simplemente los ha dejado desorientados. En este punto, no estoy seguro de si confiaría incluso en declaraciones documentadas sobre su enfoque actual: sería cauteloso. Cuesta poco comprobar la orientación después de haber procesado todos los bordes en un anillo.
whuber
... bueno, afortunadamente puse ese pequeño descargo de responsabilidad del IIRC. ;) Ya no es tan frecuente que haga cosas de nivel de geometría en la pila de ESRI.
Dan S.
@Dan S. - Todavía no estoy claro acerca de cómo se puede programar esto en Python. ¿Hay alguna guía para esto?
djq
1

¿Puede esto ayudar?

julien
fuente
Ese es un papel interesante que desenterraste. Sin embargo, ninguno de los métodos descritos allí es relevante para un cálculo de exposición solar (a menos que el cálculo sea una aproximación cruda, lo que no parece ser el caso aquí).
whuber
1

/ * Quizás esto ayude:

Azimut - pi / 2 es la orientación hacia afuera de los lados de un polígono RHR:

Aquí hay un ejemplo de PostGIS, puede crear la tabla bldg117862 usando la instrucción al final. El SRID es EPSG 2271 (PA StatePlane North Feet) y la geometría es un Multipolígono. Para visualizar en ArcGIS 10, pegue las consultas / subconsultas en una conexión de Capa de consulta a postgis después de crear la tabla bldg117862. * /

- === INICIO DE CONSULTA ===

/ * La consulta externa proporciona orientación de las ortogonales externas y crea líneas ortogonales externas de igual longitud que las de los lados desde el punto medio de los lados.

Las direcciones de orientación dominantes serán la suma de la longitud, agrupadas por orientación, en orden descendente * /

SELECCIONE line_id como side_id, longitud, grados (ortohoz) como orientación, 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

--siguiente subconsulta externa hace líneas a partir de pares de puntos de los lados, calcula el acimut (ortoaz) de ortogonal hacia afuera para cada segmento

(SELECT bldg2009gid, line_id, st_length (st_makeline (startpoint, endpoint)) :: numeric (10,2) como length, azimuth (startpoint, endpoint), azimuth (startpoint, endpoint) - pi () / 2 as orthoaz, st_makeline ( punto de inicio, punto final) como geom de

/ * subconsulta más interna: use generate_series () para descomponer polígonos de construcción en pares de puntos de punto de inicio / punto final de los lados - nota1 - forzar la regla de la mano derecha para garantizar la orientación común de todos los lados del polígono nota2 - el ejemplo usa multipolígono, para el polígono geometryn () se puede quitar */

(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 punto de inicio, pointn (exteriorring (geometryn (st_forceRHR (geom), 1)), generate_series (2, npoints (exteriorring (geometryn (st_forceRHR (geom ), 1))))) como punto final de bldg117862) como t1) como t2

- === FIN DE LA CONSULTA ===

- la tabla bldg117862 crea / inserta instrucciones

ESTABLECER STANDARD_CONFORMING_STRINGS EN ON; SELECCIONE DropGeometryColumn ('', 'bldg117862', 'geom'); TABLA DE GOTA "bldg117862"; EMPEZAR; CREATE TABLE "bldg117862" (gid serial PRIMARY KEY, "motherpin" varchar (14), "taxpin" varchar (14), "status" varchar (15), "area" numeric, "prev_area" numeric, "pct_change" numeric, "picture" varchar (133), "mappage" varchar (6), "sref_gid" int4, "e_address" varchar (19), "a_address" varchar (19), "perim" numérico, "card" int4, "a_addnum" int4, "e_street" varchar (50), "a_street" varchar (50), "e_hsnum" varchar (10)); SELECCIONE AddGeometryColumn ('', 'bldg117862', 'geom', '2271', 'MULTIPOLYGON', 2); 0106000020DF080000010000000103000020DF080000010000000B0000008C721D6C98AC34415E2C5BB9D3E32541AE56DE17BEAC34410613E5A0A0E325411AB6C794AEAC3441BA392FE372E32541C89C38429DAC3441643857628AE325418C299A9095AC3441F66C29B573E32541983F02087EAC34413080AA9F93E325419BAC3C0A86AC3441AC1F3B3DABE32541803A40B974AC3441E8CF3DB9C2E325413E3758C186AC3441D0AAB0E7F7E325410AAAA5429BAC3441BA971217DCE325418C721D6C98AC34415E2C5BB9D3E32541' ); CREAR ÍNDICE "bldg117862_geom_gist" EN "bldg117862" utilizando gist ("geom" gist_geometry_ops); FINAL;


fuente
0

Suponiendo que la orientación de los segmentos de línea es constante en un polígono, puede calcular la demora (rumbo) de un vector perpendicular a cada segmento de línea. No tengo tiempo en este momento para descifrar el código, pero si necesita las matemáticas, se puede suministrar fácilmente :-)

WolfOdrade
fuente