Intersezione di area in Python

Question

Ho un codice che accetta una condizione C come input e calcola la soluzione al mio problema come 'area consentita' A sullo spazio (x, y). Quest'area è composta da diversi "tubi", definiti da 2 linee che non possono mai attraversare.

Il risultato finale che sto cercando deve soddisfare k condizioni {C1, .., Ck}, ed è quindi un'intersezione S tra k aree {A1, .., Ak}.

Ecco un esempio con 2 condizioni (A1: verde, 3 tubi A2: viola, 1 tubo); la soluzione S è in rosso.

Come posso trovare S quando ho a che fare con 4 aree di circa 10 tubi di ciascuno? (La trama finale è orribile!)

Avrei bisogno di essere in grado di tracciarlo, e di trovare la coordinata media e la varianza dei punti in S (varianza di ogni coordinata). [Se esiste un modo efficace per sapere se un punto P appartiene a S oppure no, userò solo un metodo Monte Carlo].

Idealmente, mi piacerebbe anche essere in grado di implementare "tubi proibiti" che rimuoverei da S [potrebbe essere un po 'più complicato di intersecare S con l'esterno della mia area proibita, poiché due tubi dal stessa area può attraversare (anche se le linee che definiscono un tubo non si incrociano mai)].

---

Nota:

• Il codice memorizza anche la lunghezza dell'arco delle linee.

• Le righe vengono memorizzate come matrici di punti (circa 1000 punti per linea). Le due linee che definiscono un tubo non hanno necessariamente lo stesso numero di punti, ma Python può interpolarle TUTTE come funzione della loro lunghezza d'arco in 1 secondo.

• Le linee sono funzioni parametriche (ad esempio, non è possibile scrivere y = f (x), poiché le linee possono essere verticali).

• La trama è stata modificata con la vernice per ottenere il risultato sulla destra ... Non molto efficiente!

---

Edit:

• non so come posso usare plt.fill_between per un incrocio multipla (posso farlo qui per 2 condizioni, ma ho bisogno del codice per farlo automaticamente quando ci sono troppe linee per il giudizio dell'occhio).

• Per ora ho appena generato le linee. Non ho scritto nulla per trovare la soluzione finale dato che non so assolutamente quale sia la struttura più adatta per questo. [Tuttavia, una versione precedente del codice è stata in grado di trovare i punti di intersezione tra le linee di 2 tubi diversi, e stavo progettando di passarli come poligoni a formosa, ma questo implicava molti altri problemi ..]

• Non penso di poterlo fare con sets: la scansione dell'intera area (x, y) alla precisione richiesta rappresenta circa 6e8 punti ...[Le linee hanno solo punti 1e3 grazie ad una dimensione di passo variabile (si adatta alla curvatura), ma l'intero problema è abbastanza grande]

Answer 1

Problema risolto con ben fatto!

Ho definito ogni tubo come Polygon e un'area A è un oggetto MultiPolygon creato come unione dei suoi tubi.

Il metodo intersection calcola quindi la soluzione che stavo cercando (la sovrapposizione tra tutte le aree).

Il tutto è quasi istantaneo. Non sapevo che fosse così bello con oggetti di grandi dimensioni [circa 2000 punti per tubo, 10 tubi per area, 4 aree].

Grazie per il vostro aiuto! :)

Modifica:

Un esempio di lavoro.

import matplotlib.pyplot as plt 
import shapely 
from shapely.geometry import Polygon 
from descartes import PolygonPatch 
import numpy as np 

def create_tube(a,height): 
    x_tube_up = np.linspace(-4,4,300) 
    y_tube_up = a*x_tube_up**2 + height 
    x_tube_down = np.flipud(x_tube_up)   #flip for correct definition of polygon 
    y_tube_down = np.flipud(y_tube_up - 2) 

    points_x = list(x_tube_up) + list(x_tube_down) 
    points_y = list(y_tube_up) + list(y_tube_down) 

    return Polygon([(points_x[i], points_y[i]) for i in range(600)]) 

def plot_coords(ax, ob): 
    x, y = ob.xy 
    ax.plot(x, y, '+', color='grey') 


area_1 = Polygon()   #First area, a MultiPolygon object 
for h in [-5, 0, 5]: 
    area_1 = area_1.union(create_tube(2, h)) 

area_2 = Polygon() 
for h in [8, 13, 18]: 
    area_2 = area_2.union(create_tube(-1, h)) 

solution = area_1.intersection(area_2)  #What I was looking for 

########## PLOT ########## 

fig = plt.figure() 
ax = fig.add_subplot(111) 

for tube in area_1: 
    plot_coords(ax, tube.exterior) 
    patch = PolygonPatch(tube, facecolor='g', edgecolor='g', alpha=0.25) 
    ax.add_patch(patch) 

for tube in area_2: 
    plot_coords(ax, tube.exterior) 
    patch = PolygonPatch(tube, facecolor='m', edgecolor='m', alpha=0.25) 
    ax.add_patch(patch) 

for sol in solution: 
    plot_coords(ax, sol.exterior) 
    patch = PolygonPatch(sol, facecolor='r', edgecolor='r') 
    ax.add_patch(patch) 

plt.show() 

E la trama: