问题描述
我有两个地理数据框-一个带有折线，另一个带有点。直到现在，我仍在处理单个多段线和几个this点。我试图以某种方式类似于描述为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']
与先前的帖子想法类似，但遇到明显的错误：
Exception has occurred: ValueError
Can only compare identically-labeled Series objects
我刚刚开始使用Python进行编程，因此也许可以采用一种非常不同的方法来解决问题。任何帮助或建议将不胜感激。
解决方法
此问题是由于数据帧没有相同的索引值引起的。在此link中，它说明了当“ .em不同时， 比较运算符引发ValueError时如何处理。 ”。可以尝试以下操作：
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,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,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,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,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,where=y1fill > y2fill,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,tof,cl in zip(xfill,under_curve,colors):
tf1 = x in x2 and y in y2
if tf1:
ax.scatter(x,c=cl) # red or green
ax.text(x,"original",color=cl,alpha=0.65)
else:
# the (additional/densified) interpolated point created by `np.interp()`
# plot them in gray color
ax.scatter(x,c='lightgray',edgecolor='gray',linewidth=0.2)
ax.text(x,"densified",color='gray',alpha=0.75)
plt.show()
输出图：
