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:

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

Fonti dati

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:

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:

splitted_line = split(line, minimum_distance)

for line in list(splitted_line):
    print(line)
LINESTRING (438260.8093535866 4528121.6968436185, 442844.8639858717 4531930.850183684)
/tmp/ipykernel_8094/1721271714.py:3: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  for line in list(splitted_line):
/tmp/ipykernel_8094/1721271714.py:3: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in list(splitted_line):

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?

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:

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.

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.

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

first_linestring2 = substring(
    geom=line,
    start_dist=line_length*0,
    end_dist=line_length*distance
)
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.

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()
/tmp/ipykernel_8094/2593387518.py:8: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in splitted_lines:

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!

Sponsorizza

Se vuoi sponsorizza il progetto andando su Sponsor this project nella home di questo repository.

Alla prossima!