Calcul de l'orientation de chaque côté du polygone à l'aide d'ArcPy?

9

Je veux examiner l'orientation de chaque ligne dans un polygone afin de pouvoir calculer leur exposition solaire. Chaque polygone représente un bâtiment et a une hauteur associée. Pour le moment, je veux juste considérer l'orientation, et je considérerai plus tard les problèmes d'ombrage.

Une approche que je pensais était de diviser le polygone en lignes et de calculer l'orientation de chaque ligne, mais la difficulté est que je dois ensuite identifier la face externe de cette ligne. Bien que la plupart des polygones soient de simples figures à quatre côtés avec des lignes droites, il y en a un petit nombre où ce n'est pas le cas (un problème que je veux juste considérer, mais que je n'ai pas encore à résoudre).

Je connais python et je prévoyais de faire tout cela à partir d'un script.

djq
la source
1
Faites les réponses à cette question: gis.stackexchange.com/questions/1886/…
Sean
@Sean - merci, oui cela aide en général. Je ne connais pas les commandes ArcGIS spécifiques qui peuvent être utilisées pour ce faire. En outre, la question de savoir ce que la face intérieure / extérieure de la ligne dans un polygone reste toujours.
djq
Quand vous dites «orientation», vous ne l'avez pas assez bien défini pour moi --- quelle est l'orientation d'un hexagone parfait, par exemple?
Dan S.
@ Dan S. Je viens de modifier ma réponse pour inclure l'orientation de chaque ligne dans le polygone.
djq
désolé répondeurs - je n'ai pas réussi à attribuer correctement la prime hier.
djq

Réponses:

6

Si vous voulez juste une orientation majoritaire, consultez la réponse de @Mapperz ci-dessus.

Sinon, comme vous le dites, vous pouvez diviser les polygones en lignes à l'aide de l' outil Polygone sur ligne . Cela ajoute un champ FID gauche et droit, où le champ est -1 s'il n'y a pas de polygone extérieur - cela peut provoquer un peu de déblayage si vos bâtiments sont adjacents ou se chevauchent.

À partir de là, vous pouvez diviser vos lignes à chaque sommet (utilisez peut-être Diviser en lignes COGO ), puis calculer les angles sur chacune des lignes (éventuellement en mettant à jour les attributs COGO ).

En supposant que votre champ d'angle soit calculé à partir du nord, l'aspect sera correct là où le left_FID est -1, et pour obtenir l'aspect lorsque le right_FID est -1, ajoutez simplement 180 °. Ensuite, sur la base du FID d'origine, vous pouvez agréger, obtenir l'aspect majoritaire en fonction de la longueur, etc.

L'outil Polygon to line est scriptable (pour autant que je sache), les outils COGO ne le sont pas, vous devrez donc trouver quelque chose vous-même.

J'espère que cela t'aides!

om_henners
la source
2
Plug sans vergogne, et un commentaire: tout d'abord, si vous n'avez pas de licence Info, code.google.com/p/boundary-generator fait plus ou moins la même chose que l'outil Polygon to Line. Deuxièmement - cela vous donne beaucoup plus d'informations que de simplement regarder une sorte d'orientation générale du polygone global ... par exemple, vous pouvez utiliser le tableau résultant pour faire des calculs d'exposition au soleil beaucoup plus précis, en tenant compte de l'aspect de chaque mur individuel .
Dan S.
5

Trouver l'orientation

Dans un script, le polygone sera disponible sous la forme d'un ensemble d'anneaux - un anneau extérieur et zéro ou plusieurs anneaux intérieurs - chaque anneau étant représenté par un tuple de vecteurs ordonné cycliquement (v [0], v 1 , ... , v [m-1], v [m] = v [0]). Chaque vecteur donne les coordonnées d'un sommet (sans que deux sommets successifs coïncident). À partir de cela, il est simple, comme d'autres l'ont souligné, d'obtenir des vecteurs normaux (c'est-à-dire des vecteurs perpendiculaires aux directions des bords):

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

L'opération "t" fait pivoter un vecteur de 90 degrés dans le sens antihoraire :

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

Seules les directions de ces vecteurs normaux comptent, alors redimensionnez-les pour qu'elles aient une longueur unitaire: un vecteur (x, y) redimensionne à (x / s, y / s) où s = Sqrt (x ^ 2 + y ^ 2) (qui est la longueur de son bord correspondant). À partir de maintenant, supposons que cela a été fait. Écrivez les composants des vecteurs normaux unitaires résultants comme

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

Discriminer l'extérieur de l'intérieur

Comme vous le remarquez, cela laisse une ambiguïté directionnelle: devrions-nous utiliser n [i] ou -n [i]? Lequel pointe vers l'extérieur? Cette question revient à trouver le degré de la carte de Gauss . Pour le calculer, vous devez additionner les angles selon lesquels les directions normales changent lorsque vous marchez autour d'un anneau. Comme les vecteurs normaux ont une longueur unitaire, le cosinus de l'angle entre deux bords successifs est

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

(Définissez n [m] = n [0].)

Le sinus de l'angle entre deux bords successifs est

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

(Notez que ces calculs ne nécessitent jusqu'à présent que des sommes, des différences et des produits.) L'application de la fonction de tangente inverse principale (ATan2) à une telle paire (cosinus, sinus) donne l'angle q_i entre -180 et 180 degrés. La somme de ces angles pour i = 0, 1, ..., n-1 produit (jusqu'à l'erreur en virgule flottante) la courbure totale de l'anneau, qui doit être un multiple de 360 ​​degrés; pour un anneau fermé non auto-intersecté, ce sera soit +360 soit -360. Dans le premier cas, le degré est 1 et dans le second cas, le degré est -1. Les normales sont toutes orientées vers l'extérieur lorsque le degré de l'anneau extérieur est +1 et que les degrés des anneaux intérieurs sont -1. Réorientez-les anneau par anneau, si nécessaire, selon cette règle. Autrement dit, si le degré d'un anneau est l'opposé de celui qui est nécessaire, annulez toutes les normales de cet anneau. Vous pouvez maintenant procéder à vos calculs d'insolation.

whuber
la source
Merci pour une réponse très complète. Lorsque vous dites: «Dans un script, le polygone sera disponible sous la forme d'un ensemble d'anneaux», comment le polygone est-il disponible? Bien que familier avec certains scripts python, je ne sais pas comment interpréter le polygone de cette manière. Je le lis lentement et j'essaye de le comprendre, mais j'ai du mal à traduire certaines des explications en pseudocode que je peux écrire.
djq
Quelques notes sur cette réponse: les API SIG typiques présenteront toujours l'extérieur d'un polygone dans le sens antihoraire, IIRC, ce qui signifie que vous n'avez pas à faire cette distinction sophistiquée de l'extérieur de l'intérieur - utilisez simplement des rotations dans le sens horaire sur les anneaux extérieurs et dans le sens antihoraire sur les trous. Une clarification pour celenius: le bit `` ensemble d'anneaux '' est la façon dont la plupart des API, y compris le python d'ArcGIS, vous permettent de décomposer un polygone en ses segments de ligne constitutifs - un anneau est une courbe fermée d'entre eux. Puisque les polygones peuvent supporter des îles et des trous, vous pouvez avoir plusieurs anneaux ...
Dan S.
@Dan J'aimerais que ce soit le cas que les SIG maintiennent des représentations aussi cohérentes. Au fil des ans, le logiciel ESRI a fluctué dans les deux sens entre des polygones orientés négativement et positivement et ne les laissant que peu orientés. À ce stade, je ne suis pas sûr de pouvoir compter même sur des déclarations documentées concernant leur approche actuelle: je serais prudent. Il est peu coûteux de vérifier l'orientation après avoir déjà traité tous les bords d'un anneau.
whuber
.. eh bien, j'ai mis ce petit avertissement de l'IIRC, heureusement. ;) Ce n'est pas si souvent que je fais des trucs de niveau géométrique sur la pile ESRI.
Dan S.
@Dan S. - Je ne sais toujours pas comment cela peut être programmé en python. Existe-t-il des directives à ce sujet?
djq
1

Cela peut- il aider?

julien
la source
C'est un papier intéressant que vous avez déterré. Cependant, aucune des méthodes décrites ici n'est pertinente pour un calcul de l'exposition solaire (sauf si le calcul est destiné à être une approximation brute, ce qui ne semble pas être le cas ici).
whuber
1

/ * Peut-être que cela aide:

Azimut - pi / 2 est l'orientation tournée vers l'extérieur des côtés d'un polygone RHR:

Voici un exemple PostGIS, vous pouvez créer la table bldg117862 à l'aide de l'instruction à la fin. Le SRID est EPSG 2271 (PA StatePlane North Feet) et la géométrie est un multipolygone. Pour visualiser dans ArcGIS 10, collez la requête / les sous-requêtes dans une connexion Query Layer à postgis après avoir créé la table bldg117862. * /

- === DÉBUT DE LA REQUÊTE ===

/ * La requête externe fournit l'orientation des orthogonales vers l'extérieur et crée des lignes orthogonales vers l'extérieur de longueurs égales à celles des côtés à partir du milieu des côtés.

La (les) direction (s) dominante (s) seront la somme de la longueur, groupée par orientation, en ordre décroissant * /

SELECT line_id as side_id, longueur, degrés (orthoaz) as orientation, st_makeline (st_setsrid (st_line_interpolate_point (geom, .5), 2271), st_setsrid (st_makepoint (st_x (st_line_interpolate_point (geom, .5)) + (longueur * (sin (sin ( orthoaz))), st_y (st_line_interpolate_point (geom, .5)) + (length * (cos (orthoaz))))), 2271)) as geom from

- la prochaine sous-requête externe crée des lignes à partir des paires de points des côtés, calcule l'azimut (orthoaz) de l'orthogonal extérieur pour chaque segment

(SELECT bldg2009gid, line_id, st_length (st_makeline (startpoint, endpoint)) :: numeric (10,2) as length, azimuth (startpoint, endpoint), azimuth (startpoint, endpoint) - pi () / 2 as orthoaz, st_makeline ( point de départ, point final) en tant que geom de

/ * sous-requête la plus interne - utilisez generate_series () pour décomposer les polygones de construction en paires de points de début / point de fin des côtés - note1 - force la règle de droite pour assurer l'orientation commune de tous les côtés du polygone note2 - l'exemple utilise le multipolygone, pour polygoner la géométrie () Peut être enlevé */

(SELECT generate_series (1, npoints (exteriorring (geometryn (st_forceRHR (geom), 1))) - 1) as line_id, gid as bldg2009gid, pointn (exteriorring (geometryn (st_forceRHR (geom), 1)), generate_series (1, npoints (exteriorring (geometryn (st_forceRHR (geom), 1))) - 1)) comme point de départ, pointn (exteriorring (geometryn (st_forceRHR (geom), 1)), generate_series (2, npoints (exteriorring (geometryn (st_forceRHR (geom) ), 1))))) comme point de terminaison de bldg117862) comme t1) comme t2

- === FIN DE LA REQUÊTE ===

- les instructions de création / insertion de la table bldg117862

METTEZ STANDARD_CONFORMING_STRINGS SUR ON; SELECT DropGeometryColumn ('', 'bldg117862', 'geom'); TABLE DROP "bldg117862"; COMMENCER; 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" numeric, "card" int4, "a_addnum" int4, "e_street" varchar (50), "a_street" varchar (50), "e_hsnum" varchar (10)); SELECT AddGeometryColumn ('', 'bldg117862', 'geom', '2271', 'MULTIPOLYGON', 2); 0106000020DF080000010000000103000020DF080000010000000B0000008C721D6C98AC34415E2C5BB9D3E32541AE56DE17BEAC34410613E5A0A0E325411AB6C794AEAC3441BA392FE372E32541C89C38429DAC3441643857628AE325418C299A9095AC3441F66C29B573E32541983F02087EAC34413080AA9F93E325419BAC3C0A86AC3441AC1F3B3DABE32541803A40B974AC3441E8CF3DB9C2E325413E3758C186AC3441D0AAB0E7F7E325410AAAA5429BAC3441BA971217DCE325418C721D6C98AC34415E2C5BB9D3E32541' ); CRÉER L'INDEX "bldg117862_geom_gist" ON "bldg117862" en utilisant gist ("geom" gist_geometry_ops); FIN;


la source
0

En supposant que l'orientation des segments de ligne est constante dans un polygone, vous pouvez calculer le relèvement (cap) d'un vecteur perpendiculaire à chaque segment de ligne. Je n'ai pas le temps pour l'instant de me baser sur le code mais si vous avez besoin de maths, il peut facilement être fourni :-)

WolfOdrade
la source