Encontrar el punto más cercano en otro GeoDataFrame utilizando GeoPandas

Question

Tengo dos geodataframes:

import geopandas as gpd
from shapely.geometry import Point

gpd1 = gpd.GeoDataFrame([['John',1,Point(1,1)],['Smith',1,Point(2,2)],['Soap',1,Point(0,2)]],columns=['Name','ID','geometry'])

gpd2 = gpd.GeoDataFrame([['Work',Point(0,1.1)],['Shops',Point(2.5,2)],['Home',Point(1,1.1)]],columns=['Place','geometry'])

y quiero encontrar el nombre del punto más cercano en gpd2 para cada fila de g pd1 :

desired_output = 

    Name  ID     geometry  Nearest
0   John   1  POINT (1 1)     Home
1  Smith   1  POINT (2 2)    Shops
2   Soap   1  POINT (0 2)     Work

He intentado que esto funcione utilizando una función lambda:

gpd1['Nearest'] = gpd1.apply(lambda row: min_dist(row.geometry,gpd2)['Place'] , axis=1)

con

def min_dist(point, gpd2):
    geoseries = some_function(...)
    return geoseries

Answer 1

La respuesta de Gene no me sirvió. Finalmente descubrí que gpd2.geometry.unary_union daba como resultado una geometría que sólo contenía unos 30.000 de mi total de aproximadamente 150.000 puntos. Para cualquier otra persona que se encuentre con el mismo problema, aquí es cómo lo resolví:

from shapely.ops import nearest_points
from shapely.geometry import MultiPoint

gpd2_pts_list = gpd2.geometry.tolist()
gpd2_pts = MultiPoint(gpd2_pts_list)
def nearest(point, gpd2_pts, gpd2=gpd2, geom_col='geometry', src_col='Place'):
     # find the nearest point
     nearest_point = nearest_points(point, gpd2_pts)[1]
     # return the corresponding value of the src_col of the nearest point
     value = gpd2[gpd2[geom_col] == nearest_point][src_col].get_values()[0]
     return value

gpd1['Nearest'] = gpd1.apply(lambda x: nearest(x.geometry, gpd2_pts), axis=1)

Answer 2

Si está interesado en un tutorial interactivo sobre el uso de shapely nearest_points o scikit-learn BallTree conjuntos de datos más grandes (en lugar de scipy KDTree ) He encontrado Automating GIS processes por la Universidad de Helsinki, Finlandia) para ser realmente útil

Su Github (las notas del curso se actualizan cada año): https://github.com/Automating-GIS-processes/site

---

• nearest_points (2020): https://autogis-site.readthedocs.io/en/latest/notebooks/L3/04_nearest-neighbour.html

  from shapely.ops import nearest_points

  def get_nearest_values(row, other_gdf, point_column='geometry', value_column="geometry"): """Find the nearest point and return the corresponding value from specified value column."""

  # Create an union of the other GeoDataFrame's geometries:
  other_points = other_gdf["geometry"].unary_union
  
  # Find the nearest points
  nearest_geoms = nearest_points(row[point_column], other_points)
  
  # Get corresponding values from the other df
  nearest_data = other_gdf.loc[other_gdf["geometry"] == nearest_geoms[1]]
  
  nearest_value = nearest_data[value_column].values[0]
  
  return nearest_value

---

• BallTree (2020): https://autogis-site.readthedocs.io/en/latest/notebooks/L3/06_nearest-neighbor-faster.html

  from sklearn.neighbors import BallTree import numpy as np

  def get_nearest(src_points, candidates, k_neighbors=1): """Find nearest neighbors for all source points from a set of candidate points"""

  # Create tree from the candidate points
  tree = BallTree(candidates, leaf_size=15, metric='haversine')
  
  # Find closest points and distances
  distances, indices = tree.query(src_points, k=k_neighbors)
  
  # Transpose to get distances and indices into arrays
  distances = distances.transpose()
  indices = indices.transpose()
  
  # Get closest indices and distances (i.e. array at index 0)
  # note: for the second closest points, you would take index 1, etc.
  closest = indices[0]
  closest_dist = distances[0]
  
  # Return indices and distances
  return (closest, closest_dist)

  def nearest_neighbor(left_gdf, right_gdf, return_dist=False): """ For each point in left_gdf, find closest point in right GeoDataFrame and return them.

  NOTICE: Assumes that the input Points are in WGS84 projection (lat/lon).
  """
  
  left_geom_col = left_gdf.geometry.name
  right_geom_col = right_gdf.geometry.name
  
  # Ensure that index in right gdf is formed of sequential numbers
  right = right_gdf.copy().reset_index(drop=True)
  
  # Parse coordinates from points and insert them into a numpy array as RADIANS
  left_radians = np.array(left_gdf[left_geom_col].apply(lambda geom: (geom.x * np.pi / 180, geom.y * np.pi / 180)).to_list())
  right_radians = np.array(right[right_geom_col].apply(lambda geom: (geom.x * np.pi / 180, geom.y * np.pi / 180)).to_list())
  
  # Find the nearest points
  # -----------------------
  # closest ==> index in right_gdf that corresponds to the closest point
  # dist ==> distance between the nearest neighbors (in meters)
  
  closest, dist = get_nearest(src_points=left_radians, candidates=right_radians)
  
  # Return points from right GeoDataFrame that are closest to points in left GeoDataFrame
  closest_points = right.loc[closest]
  
  # Ensure that the index corresponds the one in left_gdf
  closest_points = closest_points.reset_index(drop=True)
  
  # Add distance if requested 
  if return_dist:
      # Convert to meters from radians
      earth_radius = 6371000  # meters
      closest_points['distance'] = dist * earth_radius
  
  return closest_points

Answer 3

Esta solución es extremadamente ineficiente, pero debería funcionar para cualquier tipo de geometría (incluyendo gdfs de tipo de geometría mixta). Yo sólo probaría esto si sus gdfs son pequeños (mi caso de uso era un gdf con unas 2000 filas para el que quería encontrar la característica más cercana de otro gdf con unas 15 filas, y tardó unos segundos en un típico portátil de oficina). Sin embargo, la intersección de características puede hacer que se pierda, así que está advertido. Originalmente se basó en la solución de @RedM, pero en su lugar asignará el índice de la característica en gdf2 que está más cerca de gdf1

gdf1["gdf2_idx"] = gdf1.apply(
    lambda row1: gdf2.apply(
        lambda row2: row1.geometry.distance(row2.geometry), axis="columns").idxmin(),
    axis="columns"
)

Answer 4

Nota, la respuesta de Gene ( https://gis.stackexchange.com/a/222388/157928 ) no me funcionó. Me seguía dando errores por lo siguiente:

AttributeError: 'Series' object has no attribute 'get_values'

Esto me ha funcionado:

from shapely.ops import nearest_points
# unary union of the gpd2 geomtries 
pts3 = gpd2.geometry.unary_union
def near(point, pts=pts3):
    # find the nearest point and return the corresponding Place value
    nearest = (gpd2.geometry == nearest_points(point, pts)[1])
    nearest = nearest.index[nearest == True].tolist()
    return gpd2.loc[nearest].Place.iloc[0]
gpd1['Nearest'] = gpd1.apply(lambda row: near(row.geometry), axis=1)
gpd1