Tagliare una linea con un punto

Tagliare una linea con un punto non è un’attività banalissima da fare in Python, in questo nuovo post mostrerò tre possibili strade per farlo.
point
Data di Pubblicazione

18 aprile 2022

Questa volta l’argomento trattato non nasce da una attività simile fatta con QGIS come per gli altri post precedenti ma dalla soluzione ad un problema in cui mi sono imbattuto e che mi ha portato via un po’ di tempo con non pochi grattacapi!

La richiesta era semplice, all’apparenza, ma non avevo fatto i conti con i Floating Point problems. In pratica, avendo una linea ed un punto lontano da essa dovevo trovare la distanza minima tra punto e lina e spezzare la linea sul punto di contatto tra la congiungente tra punto e linea.

Prima di iniziare

Librerie

Per raggiungere l’obiettivo userò le librerie che seguono:

Codice
from matplotlib import pyplot as plt
from shapely import wkt
from shapely.ops import substring, nearest_points, snap, split
from shapely.geometry import Point, LineString

import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter('ignore')

Fonti dati

Codice
point = wkt.loads('POINT (437845.7257845374 4529752.584581757)')
line = wkt.loads('LINESTRING (438260.8093535866 4528121.6968436185, 442844.8639858717 4531930.850183684)')

fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlabel('X coordinate', fontsize=15)
ax.set_ylabel('Y coordinate', fontsize=15)

ax.scatter(
    *point.xy,
    label='Point',
    s=100,
    facecolor='orange',
    edgecolor='black'
)
plt.plot(*line.xy, label='Line', color='blue')

plt.legend()

plt.show()

Premessa

Ho capito delle difficoltà legate ai floating point dopo questa esperienza:

Codice
minimum_distance = nearest_points(point, line)[1]

Con nearest_points() ottengo il punto sulla linea più prossimo al mio punto di riferimento, a questo punto è naturale tagliare la linea con questo punto:

Codice
splitted_line = split(line, minimum_distance)

for line in list(splitted_line):
    print(line)
LINESTRING (438260.8093535866 4528121.6968436185, 442844.8639858717 4531930.850183684)

La funzione split() restiuisce una geometry collection composta da due LineString. Almeno questo è quello che mi sarei aspettato leggendo anche la documentazione della libreria usata, ma purtroppo anzicchè ottenere due linee ne ho ottenuta una, la linea in ingresso. Quindi mi è sorto un dubbio: ma il punto con cui ho provato a tagliare la linea è realmente sulla linea?

Codice
snap_point = snap(minimum_distance, line, 1e-8)

assert minimum_distance.wkt == snap_point.wkt
assert line.distance(minimum_distance) < 1e-8

Gli assert non hanno dato esito negativo, quindi il punto è sulla linea. E allora perchè non riesco a tagliarla??

Facendo un po’ di ricerche ho capito che il problema è legato ai floating point e così mi sono dato da fare per trovare una soluzione senza però trascurare altre possibili strade nel dubbio che magari potesse esserci un’altra strada o più di una. Infatti chiedendo un po’ in giro ho trovato altre due soluzioni oltre la mia.

1. Soluzione mia

Uso nearest_points() per trovare il punto più prossimo alla linea, quindi prendo il from_node ed il to_node della mia linea e li uso per creare due linee:

Codice
first_linestring = LineString([Point(line.coords[0]), minimum_distance])
last_linestring = LineString([minimum_distance, Point(line.coords[1])])
minimum_distance_line = LineString([point, minimum_distance])

In questo modo ho tagliato la linea ed ho creato la linea congiungente il mio punto iniziale con la linea stessa.

Codice
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlabel('X coordinate', fontsize=15)
ax.set_ylabel('Y coordinate', fontsize=15)

plt.plot(*first_linestring.xy, label='Line', color='blue')
plt.plot(*last_linestring.xy, label='Splitted line', color='green')
plt.plot(*minimum_distance_line.xy, label='Minimum distance line', color='grey')
ax.scatter(
    *point.xy,
    label='Point',
    s=100,
    facecolor='orange',
    edgecolor='black'
)
ax.scatter(
    *minimum_distance.xy,
    label='Point on linestring',
    s=100,
    facecolor='red',
    edgecolor='black'
)

plt.legend()

plt.show()

2. Soluzione basata su projection()

Soluzione suggerita da due utenti sui social. In pratica si proietta il punto sulla linea andando ad individuare il punto di proiezione sulla stessa con cui si costruiscono le linee risultanti dal taglio.

Codice
distance = line.project(point, normalized=True)
distance
0.12131606906976151
Codice
line_length = line.length

first_linestring2 = substring(
    geom=line,
    start_dist=line_length*0,
    end_dist=line_length*distance
)
Codice
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlabel('X coordinate', fontsize=15)
ax.set_ylabel('Y coordinate', fontsize=15)

ax.scatter(
    *point.xy,
    label='Point',
    s=100,
    facecolor='orange',
    edgecolor='black'
)
plt.plot(*line.xy, label='Line', color='blue')
plt.plot(*first_linestring.xy, label='Splitted line', color='green')

plt.legend()

plt.show()

L’altra linea ed il punto di taglio si estraggono facilmente con procedimenti simili a quelli precedenti.

3. Soluzione basata su snap()

Soluzione suggerita da M.Laloux qui. In realtà ce ne è anche un’altra basata su buffer() ma avevo già incontrato precedetemente problemi legati ai floating point con quella funzione per cui l’ho scartata. La funzione usata per questa soluzione aggancia il punto di taglio alla linea e quindi effettua il taglio.

Codice
splitted_lines = split(snap(line, minimum_distance, 0.0001), minimum_distance)

fig, ax = plt.subplots(figsize=(10, 10))
ax.set_xlabel('X coordinate', fontsize=15)
ax.set_ylabel('Y coordinate', fontsize=15)

new_lines = [minimum_distance_line]
for line in splitted_lines:
    plt.plot(*line.xy, label='Splitted line', color='blue')
    new_lines.append(line)
plt.plot(*minimum_distance_line.xy, label='Minimum distance line', color='grey')
ax.scatter(
    *point.xy,
    label='Point',
    s=100,
    facecolor='orange',
    edgecolor='black'
)
ax.scatter(
    *minimum_distance.xy,
    label='Point on linestring',
    s=100,
    facecolor='red',
    edgecolor='black'
)

plt.legend()

plt.show()

Conclusione

Purtroppo i problemi legati ai floating points in shapely restano irrisolti da un po’ di anni. Le soluzioni proproste si pongoo nel caso più semplice: una linea composta da soli due vertex e cioè il from_node ed il to_node. Per situazioni più complesse si potrebbe iterare il discorso andando a segmentare la linea per poi ricongiungere i segmenti che fanno parte delle linee risultati dal taglio.

PS: i riquadri di warning che si vedono ogni tanto fanno parte del gioco, sono avvertimenti per prossime modifiche da fare al codice per poter essere utilizzabile con la prossima versione della libreria di turno.

Hai commenti, indicazioni o soluzioni alternative in Python? Scrivi nella sezione Discussioni.

Newsletter

Se vuoi ricevere aggiornamenti sui prossimi articoli iscriviti alla newsletter!

Telegram

Segui il canale Telegram CaffèGIS - I GIS a supporto dei processi evolutivi territoriali per ricevere aggiornamenti su iniziative ed informazioni a tema GIS.

Offrimi una birra!

Se trovi interessanti i miei articoli offrimi una birra!

Alla prossima!