Packen von Polygonen in Polygone mit ArcGIS Desktop?

25

Ich habe ein Boolesches Raster.

In den grauen Bereichen des Rasters möchte ich ein Polygon mit einer bestimmten Größe in einem zusammenhängenden Bereich einfügen.

Grundsätzlich habe ich ein unregelmäßiges Polygon und möchte ein bekanntes Polygon so oft wie möglich in die Ausdehnung des unregelmäßigen Polygons "einpassen".

Die Richtung des Polygons spielt keine Rolle und kann ein Quadrat sein. Ich möchte, dass es grafisch passt, aber wenn es nur eine Zahl an das Polygon anfügt (# das passt), würde das auch funktionieren.

Ich verwende ArcGIS Desktop 10.

Thad
quelle
8
Das ist ein sehr schweres Problem. Es erfordert viel Arbeit, um beispielsweise so viele Kreise wie möglich in ein Quadrat einzufügen. Wenn das ursprüngliche Polygon wie in der Abbildung kompliziert ist, benötigen Sie einige leistungsstarke Optimierungsverfahren. Die beste Methode, die ich für dieses Problem gefunden habe, ist das simulierte Tempern, das in ArcGIS jedoch nicht verfügbar ist, und das Skripting ist äußerst aufwendig (ArcGIS ist zu langsam). Könnten Sie vielleicht Ihre Anforderungen ein wenig lockern, z. B. das kleinere Polygon ausreichend oft anpassen , anstatt so oft wie möglich?
whuber
1
@whuber Danke, dass du meinen Beitrag bearbeitet hast. Ja, ausreichend oft würde funktionieren. Oder wie wäre es mit einer bestimmten Winkelausrichtung? Ex. In der obigen Abbildung habe ich das Polygon so oft angepasst, wie ich konnte, wenn ich es um 90 Grad gedreht hätte, könnten Sie noch eine weitere
Thad
1
Ja, aber es ist auch voller Fallstricke. Einige sind elementar. Beispielsweise enthielt der von ESRI verfasste und veröffentlichte Text "Einführung in ArcView GIS" (für Version 3) eine Übung, bei der ein Rechteck, das ein Fußballfeld darstellt, interaktiv in einem Polygon platziert wurde. Das Problem war, dass die Antwort der Übung falsch war, weil der Autor die Daten nicht projizierte und die Fehler bei der Verwendung der geografischen Koordinaten groß genug waren, um das Ergebnis zu beeinflussen. Die Antwort im GIS sah gut aus, aber wenn jemand versucht hätte, dieses Feld aufzubauen, hätte er festgestellt, dass nicht genug Platz dafür vorhanden war :-).
Whuber
6
@whuber Ich denke, sie dachten, eine "Ball Park" -Figur wäre ausreichend.
Kirk Kuykendall
2
Im allgemeinen Fall eines unregelmäßigen Polygons innerhalb eines unregelmäßigen Polygons ist dies ein rechenintensives Problem: Das Finden einer optimalen Lösung ist nicht in allen Fällen ein plausibles Ziel, und es ist aus technischer Sicht wahrscheinlich NP-vollständig: Welche Fälle dies sind, kann nicht vorherbestimmt werden. Wenn Sie das Problem erheblich einschränken, erhalten Sie mit einigen iterativen Zufallsanpassungsalgorithmen wahrscheinlich relativ hohe Zahlen. Ich habe das Gefühl, wenn dies eine Aufgabe ist, suchen sie nicht nach der richtigen Antwort, sondern nach kreativen Ansätzen.
MappingTomorrow

Antworten:

22

Es gibt viele Möglichkeiten, um dieses Problem anzugehen. Das Rasterformat der Daten legt einen rasterbasierten Ansatz nahe. Bei der Überprüfung dieser Ansätze erscheint eine Formulierung des Problems als binäres, ganzzahliges, lineares Programm vielversprechend, da sie dem Geist vieler GIS-Standortauswahlanalysen entspricht und leicht an diese angepasst werden kann.

In dieser Formulierung werden alle möglichen Positionen und Ausrichtungen der Füllpolygone aufgelistet, die ich als "Kacheln" bezeichnen werde. Mit jedem Plättchen ist ein Maß für seine "Güte" verbunden. Ziel ist es, eine Sammlung nicht überlappender Kacheln zu finden, deren Gesamtgüte so groß wie möglich ist. Hier können wir die Güte jeder Kachel als den Bereich ansehen, den sie abdeckt. (In datenreicheren und komplexeren Entscheidungsumgebungen berechnen wir die Güte möglicherweise als eine Kombination von Eigenschaften der Zellen, die in jeder Kachel enthalten sind, Eigenschaften, die sich möglicherweise auf die Sichtbarkeit, die Nähe zu anderen Dingen usw. beziehen.)

Die Einschränkungen für dieses Problem bestehen einfach darin, dass sich keine zwei Kacheln innerhalb einer Lösung überschneiden dürfen.

Dies kann auf eine für eine effiziente Berechnung förderliche Weise etwas abstrakter gestaltet werden, indem die Zellen in dem zu füllenden Polygon (der "Region") 1, 2, ..., M aufgelistet werden . Jede Kachelplatzierung kann mit einem Indikatorvektor aus Nullen und Einsen codiert werden, sodass die Einsen den Zellen entsprechen, die von der Kachel und den Nullen an anderer Stelle abgedeckt werden. In dieser Codierung können alle Informationen, die zu einer Sammlung von Kacheln benötigt werden, durch Summieren ihrer Indikatorvektoren (Komponente für Komponente, wie üblich) ermittelt werden: Die Summe ist ungleich Null, genau dort, wo mindestens eine Kachel eine Zelle bedeckt, und die Summe ist größer Überlappung von zwei oder mehr Kacheln. (Die Summe zählt effektiv den Betrag der Überlappung.)

Noch ein wenig Abstraktion: die Menge der möglichen Fliesen Platzierungen können sich aufgezählt werden, sagen wir 1, 2, ..., N . Die Auswahl eines beliebigen Satzes von Kachelplatzierungen entspricht einem Indikatorvektor, bei dem diejenigen die zu platzierenden Kacheln bezeichnen.

Hier ist eine kleine Illustration, um die Ideen zu korrigieren . Es wird mit dem Mathematica- Code geliefert, der für die Berechnungen verwendet wird, damit die Programmierschwierigkeiten (oder das Fehlen derselben) offensichtlich werden.

Zunächst stellen wir eine Region dar, die gekachelt werden soll:

region =  {{0, 0, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}, {1, 1, 1, 1}};

Abbildung 1: Region

Wenn wir die Zellen von links nach rechts nummerieren, beginnend von oben, enthält der Indikatorvektor für die Region 16 Einträge:

Flatten[region]

{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}

Verwenden wir die folgende Kachel zusammen mit allen Drehungen um ein Vielfaches von 90 Grad:

tileSet = {{{1, 1}, {1, 0}}};

Abbildung 2: Fliese

Code zum Erzeugen von Rotationen (und Reflexionen):

apply[s_List, alpha] := Reverse /@ s;
apply[s_List, beta] := Transpose[s];
apply[s_List, g_List] := Fold[apply, s, g];
group = FoldList[Append, {}, Riffle[ConstantArray[alpha, 4], beta]];
tiles = Union[Flatten[Outer[apply[#1, #2] &, tileSet, group, 1], 1]];

(Diese etwas undurchsichtige Berechnung wird in einer Antwort unter https://math.stackexchange.com/a/159159 erläutert. Sie zeigt, dass sie einfach alle möglichen Rotationen und Reflexionen einer Kachel erzeugt und dann alle doppelten Ergebnisse entfernt.)

Angenommen, wir platzieren die Kachel wie hier gezeigt:

Abbildung 3: Platzierung der Kacheln

Die Zellen 3, 6 und 7 werden bei dieser Platzierung abgedeckt. Das wird durch den Indikatorvektor bezeichnet

{0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Wenn wir diese Kachel um eine Spalte nach rechts verschieben, wäre dieser Indikatorvektor stattdessen

{0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Die Kombination aus dem Versuch, Kacheln an beiden Positionen gleichzeitig zu platzieren, wird durch die Summe dieser Indikatoren bestimmt.

{0, 0, 1, 1, 0, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Die 2 an der siebten Position zeigt diese Überlappung in einer Zelle (zweite Reihe nach unten, dritte Spalte von links). Da wir keine Überlappung wünschen, wird vorausgesetzt, dass die Summe der Vektoren in einer gültigen Lösung keine Einträge größer als 1 enthält.

Es stellt sich heraus, dass für dieses Problem 29 Kombinationen von Ausrichtung und Position für die Kacheln möglich sind. (Dies wurde mit einer einfachen Codierung gefunden, die eine erschöpfende Suche beinhaltet.) Wir können alle 29 Möglichkeiten darstellen, indem wir ihre Indikatoren als Spaltenvektoren zeichnen . (Die Verwendung von Spalten anstelle von Zeilen ist herkömmlich.) Hier ist ein Bild des resultierenden Arrays mit 16 Zeilen (eine für jede mögliche Zelle im Rechteck) und 29 Spalten:

makeAllTiles[tile_, {n_Integer, m_Integer}] := 
  With[{ m0 = Length[tile], n0 = Length[First[tile]]},
   Flatten[
    Table[ArrayPad[tile, {{i, m - m0 - i}, {j, n - n0 - j}}],  {i, 0, m - m0}, {j, 0, n - n0}], 1]];
allTiles = Flatten[ParallelMap[makeAllTiles[#, ImageDimensions[regionImage]] & , tiles], 1];
allTiles = Parallelize[
   Select[allTiles, (regionVector . Flatten[#]) >= (Plus @@ (Flatten[#])) &]];
options = Transpose[Flatten /@ allTiles];

Abbildung 4: Optionsfeld

(Die beiden vorherigen Indikatorvektoren werden in den ersten beiden Spalten links angezeigt.) Der scharfäugige Leser hat möglicherweise mehrere Möglichkeiten für die parallele Verarbeitung festgestellt: Diese Berechnungen können einige Sekunden dauern.

Das alles kann kompakt mit der Matrixnotation wiedergegeben werden:

  • F ist dieses Optionsfeld mit M Zeilen und N Spalten.

  • X ist die Anzeige eines Satzes von Fliesen Platzierungen, der Länge N .

  • b ist ein N- Vektor von Einsen.

  • R ist der Indikator für die Region; Es ist ein M- Vektor.

Die mit jeder möglichen Lösung X verbundene "Güte" entspricht RFX , da FX der Indikator für die von X abgedeckten Zellen ist und das Produkt mit R diese Werte summiert. (Wir könnten R gewichten, wenn wir möchten, dass die Lösungen bestimmte Gebiete in der Region bevorzugen oder meiden.) Dies soll maximiert werden. Weil wir es als ( RF ) schreiben können . X , es ist eine lineare Funktion von X : das ist wichtig. (Im folgenden Code centhält die Variable RF .)

Die Einschränkungen sind das

  1. Alle Elemente von X dürfen nicht negativ sein.

  2. Alle Elemente von X müssen kleiner als 1 sein (dies ist der entsprechende Eintrag in b );

  3. Alle Elemente von X müssen ganzzahlig sein.

Die Bedingungen (1) und (2) machen dies zu einem linearen Programm , während die dritte Anforderung es zu einem ganzzahligen linearen Programm macht.

Es gibt viele Pakete zum Lösen ganzzahliger linearer Programme, die genau in dieser Form ausgedrückt werden. Sie sind in der Lage, Werte von M und N in Zehntausende oder sogar Hunderttausende umzuwandeln. Das ist wahrscheinlich gut genug für einige reale Anwendungen.


Zur ersten Veranschaulichung habe ich mit dem Befehl von Mathematica 8 eine Lösung für das vorhergehende Beispiel berechnet LinearProgramming. (Dadurch wird eine lineare Zielfunktion minimiert . Durch Negieren der Zielfunktion wird die Minimierung leicht in die Maximierung umgewandelt.) In 0,011 Sekunden wurde eine Lösung (als Liste der Kacheln und ihrer Positionen) zurückgegeben:

b = ConstantArray[-1, Length[options]];
c = -Flatten[region].options;
lu = ConstantArray[{0, 1}, Length[First[options]]];
x = LinearProgramming[c, -options, b, lu, Integers, Tolerance -> 0.05];
If[! ListQ[x] || Max[options.x] > 1, x = {}];
solution = allTiles[[Select[x Range[Length[x]], # > 0 &]]];

Abbildung 5: Lösung

Die grauen Zellen befinden sich überhaupt nicht in der Region. Die weißen Zellen waren von dieser Lösung nicht bedeckt.

Sie können (von Hand) viele andere Fliesen ausarbeiten, die genauso gut sind wie diese - aber Sie können keine besseren finden. Dies ist eine mögliche Einschränkung dieses Ansatzes: Es gibt Ihnen eine beste Lösung, auch wenn es mehr als eine gibt. (Es gibt einige Problemumgehungen: Wenn wir die Spalten von X neu anordnen , bleibt das Problem unverändert, aber die Software wählt häufig eine andere Lösung aus. Dieses Verhalten ist jedoch nicht vorhersehbar.)

Als zweites Beispiel betrachten wir die Region in der Frage, um realistischer zu sein. Durch den Import und das Resampling des Bildes habe ich es mit einem Raster von 69 mal 81 dargestellt:

Abbildung 6: Region

Die Region umfasst 2156 Zellen dieses Gitters.

Um die Dinge interessant zu machen und die Allgemeinheit des linearen Programmieraufbaus zu veranschaulichen, versuchen wir, einen Großteil dieses Bereichs mit zwei Arten von Rechtecken abzudecken :

Abbildung 7: Fliesen

Einer ist 17 mal 9 (153 Zellen) und der andere ist 15 mal 11 (165 Zellen). Wir bevorzugen vielleicht die zweite, weil sie größer ist, aber die erste ist dünner und passt an engere Stellen. Wir werden sehen!

Das Programm beinhaltet jetzt N = 5589 mögliche Platzierungen von Kacheln. Es ist ziemlich groß! Nach 6,3 Sekunden Berechnungszeit hat Mathematica diese Zehn-Kacheln-Lösung entwickelt:

Abbildung 8: Lösung

Aufgrund des Durchhangs (z. B. könnten wir die untere linke Kachel um bis zu vier Spalten nach links verschieben) gibt es offensichtlich einige andere Lösungen, die sich geringfügig von dieser unterscheiden.

whuber
quelle
1
Eine frühere Version dieser Lösung (jedoch nicht ganz so gut) ist auf der Mathematica- Website unter mathematica.stackexchange.com/a/6888 verfügbar . Es kann auch erwähnenswert sein, dass eine geringfügige Variation der Formulierung verwendet werden kann, um das Problem der vollständigen Bedeckung der Region mit so wenig Fliesen wie möglich zu lösen (was natürlich einige Überlappungen zulässt): Dies würde das "Schlaglochflicken" lösen. Problem.
whuber
1
Im Interesse des Weltraums werden in dieser Antwort einige potenziell hilfreiche Verbesserungen nicht beschrieben. Nachdem Sie beispielsweise alle möglichen Kachelpositionen (als Indikatorvektoren) gefunden haben, können Sie sie alle addieren, um herauszufinden, welche Zellen tatsächlich von einer Kachel bedeckt werden können. Der Satz solcher Zellen zerfällt im zweiten Beispiel in zwei separate verbundene Komponenten. Dies bedeutet, dass das Problem in den beiden Komponenten unabhängig voneinander gelöst werden kann, wodurch die Größe (und damit die Rechenzeit) erheblich verringert wird. Solche anfänglichen Vereinfachungen sind in der Regel wichtig, um Probleme der realen Welt anzugehen.
Whuber
Großer Aufwand und Antwort. Chris Antwort war auch hilfreich. Vielen Dank an alle für die Hilfe! Funktioniert und hat mich wieder in die richtige Richtung gebracht.
Thad
Wow! Ich war an einem ähnlichen Problem interessiert und dieser Beitrag gab mir eine neue Perspektive. Vielen Dank. Was ist, wenn R größer ist (z. B. 140 x 140 x 20000), gibt es Möglichkeiten, die Berechnungskosten zu senken? Kennen Sie Papiere zu diesem Problem? Meine Suchbegriffe führen mich (bis jetzt) ​​nicht in die richtige Richtung.
Nimcap
@nimcap Dies ist eine wichtige Klasse von Problemen, so viel Forschung geht weiter. Schlüsselwörter, nach denen gesucht werden soll, beginnen mit "Mixed Integer Linear Program" und verzweigen sich von dort basierend auf dem, was Sie finden.
Whuber
5

Der Link zu " Über genetische Algorithmen für das Packen von Polygonen" in meiner Antwort auf eine ähnliche Frage unter " Algorithmus suchen", um die maximale Anzahl von Punkten innerhalb eines eingeschränkten Bereichs in einem minimalen Abstand zu platzieren? könnte nützlich sein. Es scheint, als könnte die Methode verallgemeinert werden, um mit beliebigen Containerformen (und nicht nur Rechtecken) zu arbeiten.

Kirk Kuykendall
quelle
Diese Arbeit enthält einige gute Ideen (+1), aber alle Algorithmen konzentrieren sich grundsätzlich auf das Packen von Polygonen in rechteckige Bereiche. Dies ist , weil es Packungen mit einer diskreten Datenstruktur (eine Folge von Polygonen zusammen mit ihren Orientierungen) darstellt , die eine Reihe von Verfahren , in dem der Polygone repräsentiert werden geschoben , zum Quadrat der Seiten parallel, in Richtung einer bestimmte Ecke. Es scheint, dass eine solche einfache diskrete Codierung für kompliziertere Regionen weniger effektiv wäre. Vielleicht würde eine anfängliche Vereinfachung der Regionen im Netz helfen.
whuber
2

Für die von Ihnen erwähnte stark eingeschränkte Teilmenge (quadratisches / dreieckiges Kacheln in einem Schlagloch) sollte dieser Pseudocode unter der Annahme der oben genannten expliziten Optimierungen eine ungefähre Antwort finden, indem er Sie einfach durch die Möglichkeiten mit einer hohen Auflösung führt und das Problem brachial erzwingt. Dies funktioniert nicht ordnungsgemäß in Situationen, in denen bei der Drehung einzelner Kacheln Zuwächse auftreten können, z. B. Rechteckkacheln oder ein sehr unregelmäßiger Behälter. Dies ist 1 Million Iterationen, Sie können bei Bedarf mehr versuchen.

Nehmen Sie ein Quadrat mit Seiten der Länge L an

Erstellen Sie ein Schachbrettmuster aus Quadraten, das mindestens den Abmessungen der Ausdehnung des Containers und mindestens 1 l in jede Richtung entspricht.

N = 0

DX = 0

DY = 0

DR = 0

Setzen Sie die Schachbrettposition auf den ursprünglichen Schwerpunkt zurück

Für (R = 1: 100)

Für (Y = 1: 100)

Für (X = 1: 100)

M = Anzahl der Quadrate zählen, die sich vollständig im Container befinden

Wenn (M> N)

DR = R

DY = Y

DX = X

N = M

Schachbrett um L / 100 nach Osten bewegen

Setzen Sie das Schachbrett nach Osten zurück

Schachbrett um L / 100 nach Norden bewegen

Setzen Sie die Schachbrett-Nordierung zurück

Schachbrett um 3,6 Grad im Uhrzeigersinn um den Schwerpunkt drehen

DY = DY * L

DX = DX * L

Setzen Sie das Schachbrett in die ursprüngliche Position und Drehung zurück

DR & "," & DX & "und" & DY & "sind die endgültige Übersetzungs- / Rotationsmatrix"

Schachbrett um DR drehen

Übersetzen Sie das Schachbrett von DX, DY

Wähle Quadrate aus, die sich vollständig im Container befinden

Quadrate exportieren

MappingTomorrow
quelle
Wenn Sie diese Prozedur in einer 2 x 5-Region versuchen, in der eine Zelle entlang der Mitte einer langen Kante fehlt, können Sie nur ein 2 x 2-Quadrat einfügen. Es passen jedoch leicht zwei solche Quadrate. Das Problem ist, dass sie nicht Teil eines regulären "Schachbrett" -Musters sind. Diese Schwierigkeit ist eines der Dinge, die dieses Problem ziemlich schwierig machen.
whuber
1
Jep. Wenn Sie eine Behälterform haben, die unregelmäßig genug ist, um mehrere nicht zusammenhängende reguläre Muster in der Größenordnung von jeweils wenigen Zellen zu unterstützen, ist dies alles andere als optimal. Das Hinzufügen solcher Elemente zum möglichen Speicherplatz erhöht die Verarbeitungszeit sehr schnell und erfordert ein gewisses Maß an Planung für den speziellen Fall, auf den Sie abzielen.
MappingTomorrow