Packen von Polygonen in Polygone mit ArcGIS Desktop?

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.

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}};

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}}};

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:

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];

(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 &]]];

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:

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 :

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:

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.

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.

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