如何区分地理定位折线下方和上方的点？
How to distinguish points below and above geolocalised polyline?
我有两个地理数据框 - 一个带有折线，第二个包含点；到目前为止，我一直在使用单条折线和几个看起来像 this. I am trying to achieve an outcome somehow similiar to the one decribed here.
的点
import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import pandas as pd
import geopandas as gpd
from shapely import wkt
roadpd = pd.read_csv('roads.csv',delimiter=';')
roadpd['geom'] = roadpd['geom'].apply(wkt.loads)
roadgpd = gpd.GeoDataFrame(roadpd, geometry=roadpd['geom'])
print(roadgpd)
addresspd = pd.read_csv('addresses.csv',delimiter=';')
addresspd['geom'] = addresspd['geom'].apply(wkt.loads)
addressgpd = gpd.GeoDataFrame(addresspd, geometry=addresspd['geom'])
print(addressgpd)
ROA_ID geom geometry
0 23 LINESTRING (16.8978861 52.417438, 16.8971243 5... LINESTRING (16.89789 52.41744, 16.89712 52.417...
ADG_ID geom geometry
0 6641801 POINT (16.89397537268492 52.4186416491154) POINT (16.89398 52.41864)
1 6641802 POINT (16.89458531052842 52.41814952579504) POINT (16.89459 52.41815)
2 6641803 POINT (16.89499192732443 52.41803648483478) POINT (16.89499 52.41804)
3 6641804 POINT (16.89532584370729 52.41794434305021) POINT (16.89533 52.41794)
4 6641805 POINT (16.89553809175901 52.41786913837524) POINT (16.89554 52.41787)
5 6641806 POINT (16.89429797664177 52.41856009093589) POINT (16.89430 52.41856)
6 6641807 POINT (16.89397725037543 52.41832458465358) POINT (16.89398 52.41832)
7 6641808 POINT (16.89376733989525 52.41815994989702) POINT (16.89377 52.41816)
8 6641809 POINT (16.89507474142366 52.418304817806) POINT (16.89507 52.41830)
9 6641810 POINT (16.89417332654704 52.41827195239621) POINT (16.89417 52.41827)
10 6641811 POINT (16.89432223101175 52.41829509557626) POINT (16.89432 52.41830)
到目前为止我想不出任何合理的解决方案。
under_curve = addressgpd['geom'] < roadgpd['geom']
类比之前的post想法遇到一个明显的错误：
Exception has occurred: ValueError
Can only compare identically-labeled Series objects
我刚开始在 Python 编程，所以也许采用一种非常不同的方法来解决这个问题可能会有所帮助。任何帮助或建议将不胜感激。
问题是因为数据帧没有相同的索引值。在此link中，它解释了当“Comparison operators raise ValueError when .index are different.”时如何处理以进行比较操作在忽略 .index 的同时，您可以尝试以下操作：
under_curve = addressgpd['geom'].values < roadgpd['geom'].values
您发现错误的原因在其他答案中Mel已经解释过了。要修复它需要一些数据排序和致密化技术，我将在下面提供的代码中以注释的形式进行解释。
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely import wkt
#from cartopy import crs as ccrs
def xy_to_linecode(x1, y1):
"""
This creates a WKT linestring code
"""
line_code = "LINESTRING("
for ix, (x,y) in enumerate(zip(x1,y1)):
line_code += str(x) + " " + str(y)
if ix<len(x1)-1:
line_code += ", "
if ix==len(x1)-1:
line_code += ")"
return line_code
def xy_to_pntlist(x2, y2):
"""
This creates a list of WKT point codes
"""
pnt_head = "POINT("
all_pnts = []
for ix, (x,y) in enumerate(zip(x2,y2)):
pnt_code = pnt_head + str(x) + " " + str(y) + ")"
all_pnts.append(pnt_code)
return all_pnts
# Part 1 (make-up data and create geoDataFrames)
# ------
# create data for a line-string
# x1 must be sorted small to large
x1 = [15, 17, 19]
y1 = [52, 51, 54]
d1 = {'col1': [1,],
'wkt': [
xy_to_linecode(x1, y1)
#'LINESTRING (15 52, 17 51, 19 54)',
]
}
# this will be 1 row dataframe
d1f = pd.DataFrame( data=d1 )
# make a geodataframe from d1f (1 LineString)
# `shapely.wkt.loads` converts WKT code to geometry
d1f['geom'] = d1f['wkt'].apply(wkt.loads)
d1f_geo = gpd.GeoDataFrame(d1f, geometry=d1f['geom'])
# create dataframe for points
# x2 must be sorted small to large
x2 = [15.2, 16.0, 16.8, 17.2, 18.0]
y2 = [52, 51, 52, 53, 51]
ids = [101, 102, 103, 104, 105]
d2 = {'id': ids,
'wkt': xy_to_pntlist(x2, y2)
#[ 'POINT (15.2 52)',
# 'POINT (16.0 51)',
# 'POINT (16.8 52)',
# 'POINT (17.2 53)',
# 'POINT (18.0 51)']
}
d2f = pd.DataFrame( data=d2 )
# make a geodataframe from d2f (several Points)
# `shapely.wkt.loads` converts WKT code to geometry
d2f['geom'] = d2f['wkt'].apply(wkt.loads)
d2f_geo = gpd.GeoDataFrame(d2f, geometry=d2f['geom'])
# Since (x1, y1) and (x2, y2) are intentionally avoided to use
# assuming working with geodataframes from this steps on ...
# we must get (x1, y1) and (x2, y2) from the dataframes
# Part 2 (make use of the geoDataFrames)
# ------
# This will get (x1, y1) from d1f_geo, equivalent with (x1, y1) of the original data
xss = []
yss = []
for ea in d1f_geo.geometry.values:
# expect d1f_geo to have 'LineString' only
xs, ys = [], []
if ea.geom_type == 'LineString':
xs.append(ea.xy[0])
ys.append(ea.xy[1])
#print(xs)
#print(ys)
xss.append(xs[0])
yss.append(ys[0])
# intentionally avoid using original (x1, y1)
x1 = xss[0]
y1 = yss[0]
# this will get (x2, y2) from d2f_geo, equivalent with (x2, y2) from the original data
x2s, y2s = [], []
for ea in d2f_geo.geometry.values:
# expect d2f_geo to cobtain 'POINT' only
x2s.append(ea.x)
y2s.append(ea.y)
#print(ea.x, ea.y)
# intentionally avoid using original (x2, y2)
x2 = x2s
y2 = y2s
# Part 3 (visualization)
# ------
# if Cartopy is used to plot the data, ...
# crs1 = ccrs.PlateCarree()
# fig, ax = plt.subplots(subplot_kw={'projection': crs1})
# fig, ax = plt.subplots(1, 1, figsize=(8, 5), subplot_kw={'projection': crs1})
# note: axes' labels will not be plotted in this case
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
# plot the single line in `d1f_geo`
kwargs={'linewidth': 1.5, 'linestyles': '--', 'alpha': 0.6}
d1f_geo.plot(color='black', ax=ax, **kwargs)
# xfill will contain all positions in x1 and x2
# that is, x1 or x2 is subset of xfill
xfill = np.sort(np.concatenate([x1, x2]))
# densify y1, y2 data to get equal shape data
y1fill = np.interp(xfill, x1, y1) # for the lineString's vertices
y2fill = np.interp(xfill, x2, y2) # for the series of points
# at this state, shapes of the arrays are equal,
# xfill.shape: (8,)
# y1fill.shape: (8,)
# y2fill.shape: (8,)
# the densified `line` data (xfill, y1fill) now have 8 vertices
# the densified `points` data (xfill, y2fill) also have 8 points
# filled color is painted here to visualize how the points/vertices are created
# if you dont need it, change True to False
if True:
ax.fill_between(xfill, y1fill, y2fill, where=y1fill < y2fill, interpolate=True, color='dodgerblue', alpha=0.1)
ax.fill_between(xfill, y1fill, y2fill, where=y1fill > y2fill, interpolate=True, color='crimson', alpha=0.1)
# plot numbers at linestring's vertices
# if you dont need this, change True to False
if True:
id = 0
for x,y in zip(xfill, y1fill):
id += 1
ax.text(x,y,str(id), color='blue')
# point is under or above the line ?
# this only works if `y2fill` and `y1fill` have equal shape
under_curve = None
if y2fill.shape == y1fill.shape:
under_curve = y2fill < y1fill
# prep list of colors, red (above line), green (below line)
colors = np.array(['green' if b else 'red' for b in under_curve])
# select only original points from (xfill, y2fill) to plot as green and red dots
# the extra (densified) points are plotted in gray
# label: "original" or "densified" are plotted at the points
for x,y,tof,cl in zip(xfill, y2fill, under_curve, colors):
tf1 = x in x2 and y in y2
if tf1:
ax.scatter(x, y, c=cl) # red or green
ax.text(x, y, "original", color=cl, alpha=0.65)
else:
# the (additional/densified) interpolated point created by `np.interp()`
# plot them in gray color
ax.scatter(x, y, c='lightgray', edgecolor='gray', linewidth=0.2)
ax.text(x, y, "densified", color='gray', alpha=0.75)
plt.show()
输出图：
我有两个地理数据框 - 一个带有折线，第二个包含点；到目前为止，我一直在使用单条折线和几个看起来像 this. I am trying to achieve an outcome somehow similiar to the one decribed here.
的点import numpy as np
import math
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import pandas as pd
import geopandas as gpd
from shapely import wkt
roadpd = pd.read_csv('roads.csv',delimiter=';')
roadpd['geom'] = roadpd['geom'].apply(wkt.loads)
roadgpd = gpd.GeoDataFrame(roadpd, geometry=roadpd['geom'])
print(roadgpd)
addresspd = pd.read_csv('addresses.csv',delimiter=';')
addresspd['geom'] = addresspd['geom'].apply(wkt.loads)
addressgpd = gpd.GeoDataFrame(addresspd, geometry=addresspd['geom'])
print(addressgpd)
ROA_ID geom geometry
0 23 LINESTRING (16.8978861 52.417438, 16.8971243 5... LINESTRING (16.89789 52.41744, 16.89712 52.417...
ADG_ID geom geometry
0 6641801 POINT (16.89397537268492 52.4186416491154) POINT (16.89398 52.41864)
1 6641802 POINT (16.89458531052842 52.41814952579504) POINT (16.89459 52.41815)
2 6641803 POINT (16.89499192732443 52.41803648483478) POINT (16.89499 52.41804)
3 6641804 POINT (16.89532584370729 52.41794434305021) POINT (16.89533 52.41794)
4 6641805 POINT (16.89553809175901 52.41786913837524) POINT (16.89554 52.41787)
5 6641806 POINT (16.89429797664177 52.41856009093589) POINT (16.89430 52.41856)
6 6641807 POINT (16.89397725037543 52.41832458465358) POINT (16.89398 52.41832)
7 6641808 POINT (16.89376733989525 52.41815994989702) POINT (16.89377 52.41816)
8 6641809 POINT (16.89507474142366 52.418304817806) POINT (16.89507 52.41830)
9 6641810 POINT (16.89417332654704 52.41827195239621) POINT (16.89417 52.41827)
10 6641811 POINT (16.89432223101175 52.41829509557626) POINT (16.89432 52.41830)
到目前为止我想不出任何合理的解决方案。
under_curve = addressgpd['geom'] < roadgpd['geom']
类比之前的post想法遇到一个明显的错误：
Exception has occurred: ValueError
Can only compare identically-labeled Series objects
我刚开始在 Python 编程，所以也许采用一种非常不同的方法来解决这个问题可能会有所帮助。任何帮助或建议将不胜感激。
问题是因为数据帧没有相同的索引值。在此link中，它解释了当“Comparison operators raise ValueError when .index are different.”时如何处理以进行比较操作在忽略 .index 的同时，您可以尝试以下操作：
under_curve = addressgpd['geom'].values < roadgpd['geom'].values
您发现错误的原因在其他答案中Mel已经解释过了。要修复它需要一些数据排序和致密化技术，我将在下面提供的代码中以注释的形式进行解释。
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely import wkt
#from cartopy import crs as ccrs
def xy_to_linecode(x1, y1):
"""
This creates a WKT linestring code
"""
line_code = "LINESTRING("
for ix, (x,y) in enumerate(zip(x1,y1)):
line_code += str(x) + " " + str(y)
if ix<len(x1)-1:
line_code += ", "
if ix==len(x1)-1:
line_code += ")"
return line_code
def xy_to_pntlist(x2, y2):
"""
This creates a list of WKT point codes
"""
pnt_head = "POINT("
all_pnts = []
for ix, (x,y) in enumerate(zip(x2,y2)):
pnt_code = pnt_head + str(x) + " " + str(y) + ")"
all_pnts.append(pnt_code)
return all_pnts
# Part 1 (make-up data and create geoDataFrames)
# ------
# create data for a line-string
# x1 must be sorted small to large
x1 = [15, 17, 19]
y1 = [52, 51, 54]
d1 = {'col1': [1,],
'wkt': [
xy_to_linecode(x1, y1)
#'LINESTRING (15 52, 17 51, 19 54)',
]
}
# this will be 1 row dataframe
d1f = pd.DataFrame( data=d1 )
# make a geodataframe from d1f (1 LineString)
# `shapely.wkt.loads` converts WKT code to geometry
d1f['geom'] = d1f['wkt'].apply(wkt.loads)
d1f_geo = gpd.GeoDataFrame(d1f, geometry=d1f['geom'])
# create dataframe for points
# x2 must be sorted small to large
x2 = [15.2, 16.0, 16.8, 17.2, 18.0]
y2 = [52, 51, 52, 53, 51]
ids = [101, 102, 103, 104, 105]
d2 = {'id': ids,
'wkt': xy_to_pntlist(x2, y2)
#[ 'POINT (15.2 52)',
# 'POINT (16.0 51)',
# 'POINT (16.8 52)',
# 'POINT (17.2 53)',
# 'POINT (18.0 51)']
}
d2f = pd.DataFrame( data=d2 )
# make a geodataframe from d2f (several Points)
# `shapely.wkt.loads` converts WKT code to geometry
d2f['geom'] = d2f['wkt'].apply(wkt.loads)
d2f_geo = gpd.GeoDataFrame(d2f, geometry=d2f['geom'])
# Since (x1, y1) and (x2, y2) are intentionally avoided to use
# assuming working with geodataframes from this steps on ...
# we must get (x1, y1) and (x2, y2) from the dataframes
# Part 2 (make use of the geoDataFrames)
# ------
# This will get (x1, y1) from d1f_geo, equivalent with (x1, y1) of the original data
xss = []
yss = []
for ea in d1f_geo.geometry.values:
# expect d1f_geo to have 'LineString' only
xs, ys = [], []
if ea.geom_type == 'LineString':
xs.append(ea.xy[0])
ys.append(ea.xy[1])
#print(xs)
#print(ys)
xss.append(xs[0])
yss.append(ys[0])
# intentionally avoid using original (x1, y1)
x1 = xss[0]
y1 = yss[0]
# this will get (x2, y2) from d2f_geo, equivalent with (x2, y2) from the original data
x2s, y2s = [], []
for ea in d2f_geo.geometry.values:
# expect d2f_geo to cobtain 'POINT' only
x2s.append(ea.x)
y2s.append(ea.y)
#print(ea.x, ea.y)
# intentionally avoid using original (x2, y2)
x2 = x2s
y2 = y2s
# Part 3 (visualization)
# ------
# if Cartopy is used to plot the data, ...
# crs1 = ccrs.PlateCarree()
# fig, ax = plt.subplots(subplot_kw={'projection': crs1})
# fig, ax = plt.subplots(1, 1, figsize=(8, 5), subplot_kw={'projection': crs1})
# note: axes' labels will not be plotted in this case
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
# plot the single line in `d1f_geo`
kwargs={'linewidth': 1.5, 'linestyles': '--', 'alpha': 0.6}
d1f_geo.plot(color='black', ax=ax, **kwargs)
# xfill will contain all positions in x1 and x2
# that is, x1 or x2 is subset of xfill
xfill = np.sort(np.concatenate([x1, x2]))
# densify y1, y2 data to get equal shape data
y1fill = np.interp(xfill, x1, y1) # for the lineString's vertices
y2fill = np.interp(xfill, x2, y2) # for the series of points
# at this state, shapes of the arrays are equal,
# xfill.shape: (8,)
# y1fill.shape: (8,)
# y2fill.shape: (8,)
# the densified `line` data (xfill, y1fill) now have 8 vertices
# the densified `points` data (xfill, y2fill) also have 8 points
# filled color is painted here to visualize how the points/vertices are created
# if you dont need it, change True to False
if True:
ax.fill_between(xfill, y1fill, y2fill, where=y1fill < y2fill, interpolate=True, color='dodgerblue', alpha=0.1)
ax.fill_between(xfill, y1fill, y2fill, where=y1fill > y2fill, interpolate=True, color='crimson', alpha=0.1)
# plot numbers at linestring's vertices
# if you dont need this, change True to False
if True:
id = 0
for x,y in zip(xfill, y1fill):
id += 1
ax.text(x,y,str(id), color='blue')
# point is under or above the line ?
# this only works if `y2fill` and `y1fill` have equal shape
under_curve = None
if y2fill.shape == y1fill.shape:
under_curve = y2fill < y1fill
# prep list of colors, red (above line), green (below line)
colors = np.array(['green' if b else 'red' for b in under_curve])
# select only original points from (xfill, y2fill) to plot as green and red dots
# the extra (densified) points are plotted in gray
# label: "original" or "densified" are plotted at the points
for x,y,tof,cl in zip(xfill, y2fill, under_curve, colors):
tf1 = x in x2 and y in y2
if tf1:
ax.scatter(x, y, c=cl) # red or green
ax.text(x, y, "original", color=cl, alpha=0.65)
else:
# the (additional/densified) interpolated point created by `np.interp()`
# plot them in gray color
ax.scatter(x, y, c='lightgray', edgecolor='gray', linewidth=0.2)
ax.text(x, y, "densified", color='gray', alpha=0.75)
plt.show()
输出图：