Obliczanie orientacji każdej strony wielokąta za pomocą ArcPy?

9

Chcę zbadać orientację każdej linii w wielokącie, aby móc obliczyć ich ekspozycję na słońce. Każdy wielokąt reprezentuje budynek i ma przypisaną wysokość. W tej chwili chcę po prostu wziąć pod uwagę orientację, a później rozważy problemy z zacienianiem.

Jednym podejściem, które myślałem, było podzielenie wielokąta na linie i obliczenie orientacji każdej linii, ale trudność polega na tym, że muszę zidentyfikować zewnętrzną powierzchnię tej linii. Chociaż większość wielokątów jest prostymi czterobocznymi figurami z prostymi liniami, istnieje niewielka liczba, w której tak nie jest (problem, który chcę po prostu rozważyć, ale jeszcze nie muszę go rozwiązać).

Znam Pythona i planowałem zrobić to wszystko ze skryptu.

djq
źródło
1
Odpowiedzi na to pytanie pomogą: gis.stackexchange.com/questions/1886/...
Sean
@Sean - dzięki, tak, to ogólnie pomaga. Nie znam konkretnych poleceń ArcGIS, których można użyć do tego celu. Ponadto pozostaje kwestia znajomości wewnętrznej / zewnętrznej powierzchni linii w wielokącie.
djq
Kiedy mówisz „orientacja”, nie dość dobrze ją dla mnie zdefiniowałeś - na przykład, jaka jest orientacja idealnego sześciokąta?
Dan S.
@ Dan S. Właśnie zredagowałem swoją odpowiedź, aby uwzględnić orientację każdej linii w wielokącie.
djq
przepraszam, odpowiadajcie - wczoraj nie udało mi się poprawnie przypisać nagrody.
djq

Odpowiedzi:

6

Jeśli chcesz orientacji większościowej, sprawdź odpowiedź @Mapperz powyżej.

W przeciwnym razie, jak mówisz, możesz podzielić wieloboki na linie za pomocą narzędzia Wielobok do linii . Dodaje to lewe i prawe pole FID, gdzie pole ma wartość -1, jeśli nie ma zewnętrznego wielokąta - może to jednak powodować pewne szarpanie się, jeśli budynki sąsiadują ze sobą lub się pokrywają.

Stamtąd możesz podzielić linie na każdym wierzchołku (być może użyj opcji Podziel na linie COGO ), a następnie obliczyć kąty na każdej z linii (potencjalnie poprzez aktualizację atrybutów COGO ).

Zakładając, że pole kąta jest obliczane od północy, aspekt będzie poprawny tam, gdzie left_FID wynosi -1, a aby uzyskać aspekt, gdy right_FID wynosi -1, po prostu dodaj 180 °. Następnie na podstawie oryginalnego FID można agregować, uzyskać większość w oparciu o długość itp.

Narzędzie „Wielokąt do linii” jest skryptowalne (o ile mi wiadomo) narzędzia COGO nie są, więc sam musiałbyś coś wymyślić.

Mam nadzieję że to pomoże!

om_henners
źródło
2
Bezwstydna wtyczka i komentarz: Po pierwsze, jeśli nie masz licencji Info, code.google.com/p/boundary-generator działa mniej więcej tak samo jak narzędzie Wielokąt do linii. Po drugie - daje to znacznie więcej informacji niż tylko spojrzenie na jakąś ogólną orientację ogólnego wielokąta ... na przykład możesz użyć tabeli wynikowej do wykonania znacznie bardziej precyzyjnych obliczeń ekspozycji na słońce, biorąc pod uwagę aspekt każdej ściany .
Dan S.
5

Znalezienie orientacji

W skrypcie wielokąt będzie dostępny jako zestaw pierścieni - jeden pierścień zewnętrzny i zero lub więcej pierścieni wewnętrznych - z każdym pierścieniem reprezentowanym przez cyklicznie uporządkowaną krotkę wektorów (v [0], v 1 , ... , v [m-1], v [m] = v [0]). Każdy wektor podaje współrzędne wierzchołka (bez dwóch kolejnych wierzchołków pokrywających się). Z tego łatwo jest, jak zauważyli inni, uzyskać normalne wektory (tj. Wektory prostopadłe do kierunków krawędzi):

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

Operacja „t” powoduje obrót wektora o 90 stopni w lewo :

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

Liczą się tylko kierunki tych normalnych wektorów, więc przeskaluj je, aby miały długość jednostkową: wektor (x, y) przeskaluje do (x / s, y / s) gdzie s = Sqrt (x ^ 2 + y ^ 2) (które jest długością odpowiedniej krawędzi). Odtąd załóżmy, że zostało to zrobione. Napisz składniki wynikowej jednostki normalnych wektorów jako

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

Dyskryminacja na zewnątrz od wewnątrz

Jak zauważyłeś, pozostawia to dwuznaczność kierunkową: czy powinniśmy używać n [i] czy -n [i]? Który wskazuje na zewnątrz? To pytanie jest równoznaczne ze znalezieniem stopnia na mapie Gaussa . Aby to obliczyć, musisz zsumować kąty, o jakie zmieniają się normalne kierunki, gdy maszerujesz wokół pierścienia. Ponieważ normalne wektory mają długość jednostkową, cosinus kąta między dwiema kolejnymi krawędziami wynosi

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

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

Sinus kąta między dwiema kolejnymi krawędziami wynosi

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

(Należy zauważyć, że te obliczenia wymagają do tej pory jedynie sum, różnic i produktów.) Zastosowanie głównej odwrotnej funkcji stycznej (ATan2) do dowolnej takiej pary (cosinus, sinus) daje kąt q_i między -180 a 180 stopni. Zsumowanie tych kątów dla i = 0, 1, ..., n-1 daje (do błędu zmiennoprzecinkowego) całkowitą krzywiznę pierścienia, która musi być wielokrotnością 360 stopni; dla zamkniętego nie-przecinającego się pierścienia będzie to +360 lub -360. W pierwszym przypadku stopień wynosi 1, aw drugim przypadku stopień -1. Wszystkie normalne są zorientowane na zewnątrz, gdy stopień pierścienia zewnętrznego wynosi +1, a stopnie pierścieni wewnętrznych - -1. W razie potrzeby ponownie ustaw je pierścień po pierścieniu, zgodnie z tą zasadą. Oznacza to, że jeśli stopień dowolnego pierścienia jest przeciwny do tego, który jest potrzebny, zaneguj wszystkie wartości normalne dla tego pierścienia. Teraz możesz przystąpić do obliczeń nasłonecznienia.

Whuber
źródło
Dziękuję za bardzo wyczerpującą odpowiedź. Kiedy powiesz: „W skrypcie wielokąt będzie dostępny jako zestaw pierścieni”, w jaki sposób dostępny jest wielokąt? Chociaż jestem zaznajomiony z niektórymi skryptami w języku Python, nie wiem, jak interpretować wielokąt w ten sposób. Czytam to powoli i staram się to zrozumieć, ale mam trudności z przetłumaczeniem niektórych wyjaśnień na pseudokod, który mogę napisać.
djq
Kilka uwag na temat tej odpowiedzi: Typowe interfejsy API GIS zawsze będą przedstawiać zewnętrzną stronę wielokąta w kolejności przeciwnej do ruchu wskazówek zegara, IIRC, co oznacza, że ​​nie musisz robić tego fantazyjnego rozróżniania na zewnątrz od wewnątrz - wystarczy użyć obrotów zgodnie z ruchem wskazówek zegara na pierścieniach zewnętrznych i na otworach przeciwnie do ruchu wskazówek zegara. Wyjaśnienie dotyczące celeniusa: bit „zestawu pierścieni” to sposób, w jaki większość interfejsów API, w tym python ArcGIS, umożliwia rozkład wieloboku na składowe segmenty linii - pierścień jest ich zamkniętą krzywą. Ponieważ wielokąty mogą obsługiwać wyspy i dziury, możesz mieć wiele pierścieni ...
Dan S.
@ Czy chciałbym, aby GISes utrzymywały tak spójne reprezentacje. Przez lata oprogramowanie ESRI zmieniało się w wielokątach negatywnie i pozytywnie, pozostawiając je niezorientowane. W tym momencie nie jestem pewien, czy oparłbym się nawet na udokumentowanych oświadczeniach dotyczących ich obecnego podejścia: byłbym ostrożny. Sprawdzanie orientacji po przetworzeniu wszystkich krawędzi w pierścieniu kosztuje niewiele.
whuber
.. cóż, na szczęście włączyłem to małe zastrzeżenie IIRC. ;) Nie zawsze tak często robię rzeczy na poziomie geometrii na stosie ESRI.
Dan S.
@ Dan S. - Nadal jestem trochę niejasny, jak można to zaprogramować w Pythonie. Czy są na to jakieś wytyczne?
djq
1

Czy to może pomóc?

Julien
źródło
To interesujący papier, który wykopałeś. Jednak żadna z opisanych tam metod nie ma zastosowania do obliczania ekspozycji na słońce (chyba że obliczenia te mają być przybliżonym przybliżeniem, co nie wydaje się tak w tym przypadku).
whuber
1

/ * Może to pomaga:

Azymut - pi / 2 to skierowana na zewnątrz orientacja boków wielokąta RHR:

Oto przykład PostGIS, możesz utworzyć tabelę bldg117862, używając instrukcji na końcu. SRID to EPSG 2271 (PA StatePlane North Feet), a geometria jest wieloboczna. Aby wizualizować w ArcGIS 10, wklej zapytanie / podkwerendy do połączenia warstwy zapytań do postgis po utworzeniu tabeli bldg117862. * /

- === START ZAPYTAŃ ===

/ * Zapytanie zewnętrzne zapewnia orientację zewnętrznych ortogonali i tworzy zewnętrzne linie ortogonalne o równej długości jak linie boków od punktu środkowego boków.

Dominujący kierunek (kierunki) będzie sumą długości, pogrupowaną według orientacji, w porządku malejącym * /

WYBIERZ identyfikator linii jako identyfikator boku, długość, stopnie (ortoaz) jako orientację, st_makeline (st_setsrid (st_line_interpolate_point (geom, .5), 2271), st_setsrid (st_makepoint (st_x (st_line_interpolate_point (geom, .5)) + + długość * (sin ( orthoaz))), st_y (point_interpolate_point (geom, .5)) + (length * (cos (orthoaz)))), 2271)) as geom from

- następne zewnętrzne podzapytanie tworzy linie z par punktowych boków, oblicza azymut (ortoaz) na zewnątrz prostopadłego dla każdego segmentu

(SELECT bldg2009gid, line_id, st_length (st_makeline (punkt początkowy, punkt końcowy)) :: numeryczne (10,2) jako długość, azymut (punkt początkowy, punkt końcowy), azymut (punkt początkowy, punkt końcowy) - pi () / 2 jako orthoaz, st_makeline ( punkt początkowy, punkt końcowy) jako geom z

/ * najbardziej wewnętrzne podzapytanie - użyj funkcji generuj_series (), aby rozłożyć budujące wielokąty na pary punkt początkowy / punkt końcowy boków - uwaga 1 - wymusza regułę prawej ręki, aby zapewnić wspólną orientację wszystkich boków wielokąta uwaga 2 - przykład używa wieloboku, dla wieloboku geometryn () może być usunięty */

(WYBIERZ wygeneruj_series (1, npoints (externalring (geometryn (st_forceRHR (geom), 1))) - 1) jako line_id, gid jako bldg2009gid, pointn (exteriorring (geometryn (st_forceRHR (geom), 1)), wygeneruj serie (1, npoints (exteriorring (geometryn (st_forceRHR (geom), 1))) - 1)) jako punkt początkowy, pointn (exteriorring (geometryn (st_forceRHR (geom), 1)), generuj serie (2, npoints (externalring (geometryn (st_forceRHR (geom) ), 1))))) jako punkt końcowy od bldg117862) jako t1) jako t2

- === KONIEC ZAPYTANIA ===

- instrukcje bldg117862 tworzą / wstawiają tabelę

ZESTAW STANDARD_CONFORMING_STRINGS NA WŁ; SELECT DropGeometryColumn („”, „bldg117862”, „geom”); TABELA UPADKU „bldg117862”; ZACZYNAĆ; UTWÓRZ TABELĘ „bldg117862” (klucz seryjny gid, „motherpin” varchar (14), „taxpin” varchar (14), „status” varchar (15), „area” numeric, „prev_area” numeric, „pct_change” numeric, „obraz” varchar (133), „mapa strony” varchar (6), „sref_gid” int4, „e_address” varchar (19), „a_address” varchar (19), „perim” numeryczny, „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' ); UTWÓRZ INDEKS „bldg117862_geom_gist„ ON ”bldg117862” przy użyciu gist („geom” gist_geometry_ops); KONIEC;


źródło
0

Zakładając, że orientacja segmentów linii jest stała w wielokącie, można obliczyć namiar (nagłówek) wektora prostopadłego do każdego segmentu linii. Nie mam teraz czasu, aby wypróbować kod, ale jeśli potrzebujesz matematyki, można ją łatwo podać :-)

WolfOdrade
źródło