-
Notifications
You must be signed in to change notification settings - Fork 2
Open
Description
| def transect_timeseries(shorelines_path, |
Dan found a bug with the transect_timeseries function that causes the multipoints not be replaced for the same shoreline that intersected multiple transects. I posted the new code that will fix this issue
Short Version
- After intersecting the shoreline with the transect the row index was duplicated
- The code to assign a replace a multipoint intersection point with a single intersection point was updating multiple rows if the shoreline row index was duplicated
New Code
# get points, keep highest cross distance point if multipoint (most seaward intersection)
joined_gdf['intersection_point'] = joined_gdf.geometry.intersection(joined_gdf['geometry_saved'])
# reset the index because running the intersection function changes the index
joined_gdf = joined_gdf.reset_index(drop=True)
for i in range(len(joined_gdf['intersection_point'])):
point = joined_gdf['intersection_point'].iloc[i]
start_x = joined_gdf['x_start'].iloc[i]
start_y = joined_gdf['y_start'].iloc[i]
if type(point) == shapely.MultiPoint:
points = [shapely.Point(coord) for coord in point.geoms]
points = gpd.GeoSeries(points, crs=crs)
coords = points.get_coordinates()
dists = [None]*len(coords)
for j in range(len(coords)):
dists[j] = cross_distance(start_x, start_y, coords['x'].iloc[j], coords['y'].iloc[j])
max_dist_idx = np.argmax(dists)
last_point = points[max_dist_idx]
# joined_gdf['intersection_point'].iloc[i] = last_point # replace this line
# This new line updates the shoreline at index (i) for a single value. We only want to update a single shoreline
joined_gdf.at[i, 'intersection_point'] = last_point # is only edits a single intersection point
Full Explanation
After intersecting the shorelines with the transects the row index is changed to the number of shorelines (0-255) and it is duplicated for every date. This means that for each instance the shoreline intersected a transect it is duplicated with the same index but a different intersection point. So it the shoreline intersected 5 transects there are 5 rows with the index 0.
joined_gdf['intersection_point'] = joined_gdf.geometry.intersection(joined_gdf['geometry_saved'])
# reset the index because running the intersection function changes the index
joined_gdf = joined_gdf.reset_index(drop=True)
The second issue was in the loop.
for i in range(len(joined_gdf['intersection_point'])):
point = joined_gdf['intersection_point'].iloc[i]
start_x = joined_gdf['x_start'].iloc[i]
start_y = joined_gdf['y_start'].iloc[i]
if type(point) == shapely.MultiPoint:
points = [shapely.Point(coord) for coord in point.geoms]
points = gpd.GeoSeries(points, crs=crs)
coords = points.get_coordinates()
dists = [None]*len(coords)
for j in range(len(coords)):
dists[j] = cross_distance(start_x, start_y, coords['x'].iloc[j], coords['y'].iloc[j])
max_dist_idx = np.argmax(dists)
last_point = points[max_dist_idx]
# This line causes all the intersection points for the index i to be changed to the last point (BAD if the index is duplicated)
# joined_gdf['intersection_point'].iloc[i] = last_point
# This new line updates the shoreline at index (i) for a single value. We only want to update a single shoreline
joined_gdf.at[i, 'intersection_point'] = last_point # is only edits a single intersection point
Metadata
Metadata
Assignees
Labels
No labels