在三个维度 returns 点 Z 中均匀插值但结果无效
shapely interpolate in three dimensions returns Point Z but invalid results
我正在尝试沿三维线使用插值。但是，.interpolate 不会考虑 Z 轴的任何变化。
LineString([(0, 0, 0), (0, 0, 1), (0, 0, 2)]).interpolate(1, normalized=True).wkt
'POINT Z (0 0 0)'
对
LineString([(0, 0, 0), (0, 1, 0), (0, 2, 0)]).interpolate(1, normalized=True).wkt
'POINT Z (0 2 0)'
我阅读了文档，它对 3D 线没有任何说明，或者在比插值文档更高级别记录了限制。
这是一个错误吗？我不敢相信我是第一个尝试这个的人。
假设没有直接的方法来完成这个，有什么关于我自己插值的建议吗？
这确实看起来像是来自 shapely 的错误。我稍微查看了源代码，我敢打赌这是 PyGEOS 的上游问题。
无论如何，这是我整理的一个小实现：
import numpy as np
import shapely
import geopandas as gpd # Only necessary for the examples, not the actual function
def my_interpolate(input_line, input_dist, normalized=False):
'''
Function that interpolates the coordinates of a shapely LineString.
Note: If you use this function on a MultiLineString geometry, it will
"flatten" the geometry and consider all the points in it to be
consecutively connected. For example, consider the following shape:
MultiLineString(((0,0),(0,2)),((0,4),(0,6)))
In this case, this function will assume that there is no gap between
(0,2) and (0,4). Instead, the function will assume that these points
all connected. Explicitly, the MultiLineString above will be
interpreted instead as the following shape:
LineString((0,0),(0,2),(0,4),(0,6))
Parameters
----------
input_line : shapely.geometry.Linestring or shapely.geometry.MultiLineString
(Multi)LineString whose coordinates you want to interpolate
input_dist : float
Distance used to calculate the interpolation point
normalized : boolean
Flag that indicates whether or not the `input_dist` argument should be
interpreted as being an absolute number or a percentage that is
relative to the total distance or not.
When this flag is set to "False", the `input_dist` argument is assumed
to be an actual absolute distance from the starting point of the
geometry. When this flag is set to "True", the `input_dist` argument
is assumed to represent the relative distance with respect to the
geometry's full distance.
The default is False.
Returns
-------
shapely.geometry.Point
The shapely geometry of the interpolated Point.
'''
# Making sure the entry value is a LineString or MultiLineString
if ((input_line.type.lower() != 'linestring') and
(input_line.type.lower() != 'multilinestring')):
return None
# Extracting the coordinates from the geometry
if input_line.type.lower()[:len('multi')] == 'multi':
# In case it's a multilinestring, this step "flattens" the points
coords = [item for sub_list in [list(this_geom.coords) for
this_geom in input_line.geoms]
for item in sub_list]
else:
coords = [tuple(coord) for coord in list(input_line.coords)]
# Transforming the list of coordinates into a numpy array for
# ease of manipulation
coords = np.array(coords)
# Calculating the distances between points
dists = ((coords[:-1] - coords[1:])**2).sum(axis=1)**0.5
# Calculating the cumulative distances
dists_cum = np.append(0,dists.cumsum())
# Finding the total distance
dist_total = dists_cum[-1]
# Finding appropriate use of the `input_dist` value
if normalized == False:
input_dist_abs = input_dist
input_dist_rel = input_dist / dist_total
else:
input_dist_abs = input_dist * dist_total
input_dist_rel = input_dist
# Taking care of some edge cases
if ((input_dist_rel < 0) or
(input_dist_rel > 1) or
(input_dist_abs < 0) or
(input_dist_abs > dist_total)):
return None
elif ((input_dist_rel == 0) or (input_dist_abs == 0)):
return shapely.geometry.Point(coords[0])
elif ((input_dist_rel == 1) or (input_dist_abs == dist_total)):
return shapely.geometry.Point(coords[-1])
# Finding which point is immediately before and after the input distance
pt_before_idx = np.arange(dists_cum.shape[0])[(dists_cum <= input_dist_abs)].max()
pt_after_idx = np.arange(dists_cum.shape[0])[(dists_cum >= input_dist_abs)].min()
pt_before = coords[pt_before_idx]
pt_after = coords[pt_after_idx]
seg_full_dist = dists[pt_before_idx]
dist_left = input_dist_abs - dists_cum[pt_before_idx]
# Calculating the interpolated coordinates
interpolated_coords = ((dist_left / seg_full_dist) * (pt_after - pt_before)) + pt_before
# Creating a shapely geometry
interpolated_point = shapely.geometry.Point(interpolated_coords)
return interpolated_point
上述函数可用于Shapely (Multi)LineStrings。下面是它应用于简单 LineString 的示例。
input_line = shapely.geometry.LineString([(0, 0, 0),
(1, 2, 3),
(4, 5, 6)])
interpolated_point = my_interpolate(input_line, 2.5, normalized=False)
print(interpolated_point.wkt)
> POINT Z (0.6681531047810609 1.336306209562122 2.004459314343183)
下面是一个使用apply方法对LineStrings的整个GeoDataFrame进行插值的例子：
line_df = gpd.GeoDataFrame({'id':[1,
2,
3],
'geometry':[input_line,
input_line,
input_line],
'interpolate_dist':[0.5,
2.5,
6.5],
'interpolate_dist_normalized':[True,
False,
False]})
interpolated_points = line_df.apply(
lambda row: my_interpolate(input_line=row['geometry'],
input_dist=row['interpolate_dist'],
normalized=row['interpolate_dist_normalized']),
axis=1)
print(interpolated_points.apply(lambda point: point.wkt))
> 0 POINT Z (1.419876550265357 2.419876550265356 3...
> 1 POINT Z (0.6681531047810609 1.336306209562122 ...
> 2 POINT Z (2.592529850263281 3.592529850263281 4...
> dtype: object
重要提示
极端情况和错误处理
请注意，我开发的函数不能很好地处理错误。在很多情况下，它只是默默地 returns 一个 None 对象。根据您的用例，您可能想要调整该行为。
多线串
上面的函数可以用在MultiLineStrings上，但是做了一些简化和假设。如果您在 MultiLineString 几何体上使用此函数，它将“展平”几何体并认为其中的所有点都是连续连接的。例如，考虑以下形状：
MultiLineString(((0,0),(0,2)),((0,4),(0,6)))
在这种情况下，函数将假定 (0,2) 和 (0,4) 之间没有间隙。相反，该函数将假定这些点都已连接。明确地说，上面的 MultiLineString 将被解释为以下形状：
LineString((0,0),(0,2),(0,4),(0,6))
有人问我“你能沿着每个轴插值而不是三个一起插值吗？”我认为答案是肯定的，这是我使用的方法。
# Upsample to 1S intervals rather than our desired interval because resample throws
# out rows that do not fall on the desired interval, including the rows we want to keep.
int_df = df.resample('1S', origin='start').asfreq()
# For each axis, interpolate to fill in NAN values.
int_df['Latitude'] = int_df['Latitude'].interpolate(method='polynomial', order=order)
int_df['Longitude'] = int_df['Longitude'].interpolate(method='polynomial', order=order)
int_df['AGL'] = int_df['AGL'].interpolate(method='polynomial', order=order)
# Now downsample to our desired frequency
int_df = int_df.resample('5S', origin='start').asfreq()
我最初以 5S 的间隔重新采样，但这导致不在间隔边界上的任何现有点被丢弃，以支持在间隔边界上的新点。对于我的用例，这很重要。如果你想要固定的时间间隔，那么你不需要先上采样再下采样。
之后，只需对三个轴分别进行插值即可。
所以，如果我开始：
我现在有：
回答形状操作函数为什么不能在 3D/Z 上运行的问题：
来自 shapely docs。 （当版本 1.8.X 是当前版本时写这个）
A third z coordinate value may be used when constructing instances,
but has no effect on geometric analysis. All operations are performed
in the x-y plane.
为了我的目的，我也需要 Z。所以正在搜索此信息以查看是否可以使用 geopandas（使用 shaply），而不是 osgeo.ogr.
我正在尝试沿三维线使用插值。但是，.interpolate 不会考虑 Z 轴的任何变化。
LineString([(0, 0, 0), (0, 0, 1), (0, 0, 2)]).interpolate(1, normalized=True).wkt
'POINT Z (0 0 0)'
对
LineString([(0, 0, 0), (0, 1, 0), (0, 2, 0)]).interpolate(1, normalized=True).wkt
'POINT Z (0 2 0)'
我阅读了文档，它对 3D 线没有任何说明，或者在比插值文档更高级别记录了限制。
这是一个错误吗？我不敢相信我是第一个尝试这个的人。
假设没有直接的方法来完成这个，有什么关于我自己插值的建议吗？
这确实看起来像是来自 shapely 的错误。我稍微查看了源代码，我敢打赌这是 PyGEOS 的上游问题。
无论如何，这是我整理的一个小实现：
import numpy as np
import shapely
import geopandas as gpd # Only necessary for the examples, not the actual function
def my_interpolate(input_line, input_dist, normalized=False):
'''
Function that interpolates the coordinates of a shapely LineString.
Note: If you use this function on a MultiLineString geometry, it will
"flatten" the geometry and consider all the points in it to be
consecutively connected. For example, consider the following shape:
MultiLineString(((0,0),(0,2)),((0,4),(0,6)))
In this case, this function will assume that there is no gap between
(0,2) and (0,4). Instead, the function will assume that these points
all connected. Explicitly, the MultiLineString above will be
interpreted instead as the following shape:
LineString((0,0),(0,2),(0,4),(0,6))
Parameters
----------
input_line : shapely.geometry.Linestring or shapely.geometry.MultiLineString
(Multi)LineString whose coordinates you want to interpolate
input_dist : float
Distance used to calculate the interpolation point
normalized : boolean
Flag that indicates whether or not the `input_dist` argument should be
interpreted as being an absolute number or a percentage that is
relative to the total distance or not.
When this flag is set to "False", the `input_dist` argument is assumed
to be an actual absolute distance from the starting point of the
geometry. When this flag is set to "True", the `input_dist` argument
is assumed to represent the relative distance with respect to the
geometry's full distance.
The default is False.
Returns
-------
shapely.geometry.Point
The shapely geometry of the interpolated Point.
'''
# Making sure the entry value is a LineString or MultiLineString
if ((input_line.type.lower() != 'linestring') and
(input_line.type.lower() != 'multilinestring')):
return None
# Extracting the coordinates from the geometry
if input_line.type.lower()[:len('multi')] == 'multi':
# In case it's a multilinestring, this step "flattens" the points
coords = [item for sub_list in [list(this_geom.coords) for
this_geom in input_line.geoms]
for item in sub_list]
else:
coords = [tuple(coord) for coord in list(input_line.coords)]
# Transforming the list of coordinates into a numpy array for
# ease of manipulation
coords = np.array(coords)
# Calculating the distances between points
dists = ((coords[:-1] - coords[1:])**2).sum(axis=1)**0.5
# Calculating the cumulative distances
dists_cum = np.append(0,dists.cumsum())
# Finding the total distance
dist_total = dists_cum[-1]
# Finding appropriate use of the `input_dist` value
if normalized == False:
input_dist_abs = input_dist
input_dist_rel = input_dist / dist_total
else:
input_dist_abs = input_dist * dist_total
input_dist_rel = input_dist
# Taking care of some edge cases
if ((input_dist_rel < 0) or
(input_dist_rel > 1) or
(input_dist_abs < 0) or
(input_dist_abs > dist_total)):
return None
elif ((input_dist_rel == 0) or (input_dist_abs == 0)):
return shapely.geometry.Point(coords[0])
elif ((input_dist_rel == 1) or (input_dist_abs == dist_total)):
return shapely.geometry.Point(coords[-1])
# Finding which point is immediately before and after the input distance
pt_before_idx = np.arange(dists_cum.shape[0])[(dists_cum <= input_dist_abs)].max()
pt_after_idx = np.arange(dists_cum.shape[0])[(dists_cum >= input_dist_abs)].min()
pt_before = coords[pt_before_idx]
pt_after = coords[pt_after_idx]
seg_full_dist = dists[pt_before_idx]
dist_left = input_dist_abs - dists_cum[pt_before_idx]
# Calculating the interpolated coordinates
interpolated_coords = ((dist_left / seg_full_dist) * (pt_after - pt_before)) + pt_before
# Creating a shapely geometry
interpolated_point = shapely.geometry.Point(interpolated_coords)
return interpolated_point
上述函数可用于Shapely (Multi)LineStrings。下面是它应用于简单 LineString 的示例。
input_line = shapely.geometry.LineString([(0, 0, 0),
(1, 2, 3),
(4, 5, 6)])
interpolated_point = my_interpolate(input_line, 2.5, normalized=False)
print(interpolated_point.wkt)
> POINT Z (0.6681531047810609 1.336306209562122 2.004459314343183)
下面是一个使用apply方法对LineStrings的整个GeoDataFrame进行插值的例子：
line_df = gpd.GeoDataFrame({'id':[1,
2,
3],
'geometry':[input_line,
input_line,
input_line],
'interpolate_dist':[0.5,
2.5,
6.5],
'interpolate_dist_normalized':[True,
False,
False]})
interpolated_points = line_df.apply(
lambda row: my_interpolate(input_line=row['geometry'],
input_dist=row['interpolate_dist'],
normalized=row['interpolate_dist_normalized']),
axis=1)
print(interpolated_points.apply(lambda point: point.wkt))
> 0 POINT Z (1.419876550265357 2.419876550265356 3...
> 1 POINT Z (0.6681531047810609 1.336306209562122 ...
> 2 POINT Z (2.592529850263281 3.592529850263281 4...
> dtype: object
重要提示
极端情况和错误处理
请注意，我开发的函数不能很好地处理错误。在很多情况下，它只是默默地 returns 一个 None 对象。根据您的用例，您可能想要调整该行为。
多线串
上面的函数可以用在MultiLineStrings上，但是做了一些简化和假设。如果您在 MultiLineString 几何体上使用此函数，它将“展平”几何体并认为其中的所有点都是连续连接的。例如，考虑以下形状：
MultiLineString(((0,0),(0,2)),((0,4),(0,6)))
在这种情况下，函数将假定 (0,2) 和 (0,4) 之间没有间隙。相反，该函数将假定这些点都已连接。明确地说，上面的 MultiLineString 将被解释为以下形状：
LineString((0,0),(0,2),(0,4),(0,6))
有人问我“你能沿着每个轴插值而不是三个一起插值吗？”我认为答案是肯定的，这是我使用的方法。
# Upsample to 1S intervals rather than our desired interval because resample throws
# out rows that do not fall on the desired interval, including the rows we want to keep.
int_df = df.resample('1S', origin='start').asfreq()
# For each axis, interpolate to fill in NAN values.
int_df['Latitude'] = int_df['Latitude'].interpolate(method='polynomial', order=order)
int_df['Longitude'] = int_df['Longitude'].interpolate(method='polynomial', order=order)
int_df['AGL'] = int_df['AGL'].interpolate(method='polynomial', order=order)
# Now downsample to our desired frequency
int_df = int_df.resample('5S', origin='start').asfreq()
我最初以 5S 的间隔重新采样，但这导致不在间隔边界上的任何现有点被丢弃，以支持在间隔边界上的新点。对于我的用例，这很重要。如果你想要固定的时间间隔，那么你不需要先上采样再下采样。
之后，只需对三个轴分别进行插值即可。
所以，如果我开始：
我现在有：
回答形状操作函数为什么不能在 3D/Z 上运行的问题：
来自 shapely docs。 （当版本 1.8.X 是当前版本时写这个）
A third z coordinate value may be used when constructing instances, but has no effect on geometric analysis. All operations are performed in the x-y plane.
为了我的目的，我也需要 Z。所以正在搜索此信息以查看是否可以使用 geopandas（使用 shaply），而不是 osgeo.ogr.