Optimisation de la distance GeoPandas

Oct 15 2020

J'ai créé la fonction suivante qui mesure la distance d'un point dans un GDF à tous les points dans un autre GDF et renvoie une table avec la distance la plus courte pour chaque point. Cela fonctionne bien pour un point mais j'ai négligé le fait d'avoir une table de 4000 points et cela prend donc 10 minutes. Je l'ai exécuté dans PostGIS et je peux le réduire en moins d'une seconde. Existe-t-il un moyen de le faire en Python qui pourrait correspondre à la vitesse de PostGIS?

def get_distance_to(gdf_in, aoi_df, aoi):
   dist_df_list = list()
   for row in range(len(gdf_in)):
      single_row = gdf_in.iloc[row]
      distances = aoi_df.geometry.distance(single_row.geometry)
      dist_list = distances.to_list()
      closest_aoi = min(dist_list)
      single_row["dist_to_"+aoi] = closest_aoi
      df = single_row.to_frame().T
      dist_df_list.append(df)
completed_distances = pd.concat(dist_df_list, ignore_index=True, sort=False)
return completed_distances

mes tables d'entrée ressemblent à ceci

et la table de sortie ressemble à ceci

Réponses

3 martinfleis Oct 15 2020 at 17:07

Pour toute opération spatiale de ce type, vous devez toujours essayer d'utiliser l'index spatial. Si vous n'êtes intéressé que par la distance minimale, ce qui suit devrait vous donner une option relativement performante.

import geopandas as gpd
from shapely.geometry import Point
import pandas as pd
import random

gdf = gpd.GeoDataFrame(geometry=[Point(random.randint(0, 1000), random.randint(0, 1000)) for _ in range(1000)])
gdf2 = gpd.GeoDataFrame(geometry=[Point(random.randint(0, 1000), random.randint(0, 1000)) for _ in range(1000)])

def get_nearest_distance(left, right, initial_buffer):
    """get distance from left to right"""
    buffered = left.buffer(initial_buffer)

    distances = []
    for i in range(len(buffered)):
        geom = buffered.geometry.iloc[i]
        query = right.sindex.query(geom)
        while query.size == 0:
            query = right.sindex.query(geom.buffer(b))
            b += initial_buffer
        distances.append(right.iloc[query].distance(left.geometry.iloc[i]).min())

    return pd.Series(distances, index=left.index)

gdf['distance_to_x'] = get_nearest_distance(gdf, gdf2, 50)

Pour 1000 à 1000 points, c'est moins d'une seconde, comparé à environ une minute de code @ gene.

Pour le rendre efficace, vous devez deviner la initial_bufferdistance à laquelle vous pensez ne sera que de quelques points. S'il n'y en a pas, il étend le tampon jusqu'à ce qu'il en atteigne.

En règle générale, si vous voulez les meilleures performances de GeoPandas, vous devez utiliser la dernière version (ce code nécessite 0.8) et les pygeos de dépendance facultatifs (https://geopandas.readthedocs.io/en/latest/getting_started/install.html#using-the-optional-pygeos-dependency), ce qui peut accélérer le code ci-dessus de l'ordre de grandeur.

2 gene Oct 15 2020 at 15:22

Itérer sur des lignes dans un (Geo) DataFrame dans (Geo) Pandas est très lent, voir Approche optimale pour itérer sur un DataFrame par exemple

L'itération dans Pandas est un anti-pattern et c'est quelque chose que vous ne devriez faire que lorsque vous avez épuisé toutes les autres options. ( Comment parcourir les lignes d'un DataFrame dans Pandas )

Vous pouvez essayer d'utiliser (Geo)DataFrame.apply()et de mettre en forme: le point le plus proche comme dans GeoPandas: trouver le point le plus proche dans une autre trame de données sans foritération (voir commentaire)

import geopandas as gpd
from shapely.geometry import Point
from shapely.ops import nearest_points
gpd1 = gpd.read_file("point1.shp") # red points
gpd2 = gpd.read_file("point2.shp") # blue points
pts3 = gpd2.geometry.unary_union
 def near(point, pts=pts3):
     # find the nearest point and return the corresponding value
     nearest = gpd2.geometry == nearest_points(point, pts)[1]
     return gpd2[nearest].id.values[0],gpd2[nearest].geometry.values[0]
 gpd1['Nearest'] = gpd1.apply(lambda row: near(row.geometry)[0], axis=1)
 gpd1['geom2'] = gpd1.apply(lambda row: near(row.geometry)[1], axis=1)
 print(gpd1)
         id         geometry         Nearest        geom2
 0   1   POINT (-0.99013 0.48096)      3     POINT (-0.77574 0.64739)
 1   2   POINT (-1.00987 0.08039)      4     POINT (-0.73060 0.10860)
 2   3   POINT (-0.71932 -0.13117)     5     POINT (-0.57827 -0.08039)
 3   4   POINT (-0.90268 -0.28914)     5     POINT (-0.57827 -0.08039)

Calculez la distance

gpd1['distance'] = gpd1.apply(lambda row: row.geometry.distance(row.geom2), axis=1)
gpd1.drop('geom2', axis=1, inplace=True)
print(gpd1)
    id         geometry             Nearest     distance
 0   1   POINT (-0.99013 0.48096)      3        0.271406
 1   2   POINT (-1.00987 0.08039)      4        0.280688
 2   3   POINT (-0.71932 -0.13117)     5        0.149905
 3   4   POINT (-0.90268 -0.28914)     5        0.385759