3

After performing some basic vector operations on an OSM street-network I have artifacts in a number of locations.

import numpy as np
import pandas as pd
import geopandas as gpd
from osgeo import gdal, ogr
import matplotlib.pyplot as plt

#get the stuff
rd = gpd.read_file('./data/roads-cput_proj.geojson')
rd.set_crs(epsg=32734, inplace=True, allow_override=True)
rd['lanes'] = rd['tags'].apply(lambda x: x.get('lanes'))
rd['lanes'] = pd.to_numeric(rd['lanes'])

#--basic cleaning
rd['bridge'] = rd['tags'].apply(lambda x: x.get('bridge'))
rd['name'] = rd['tags'].apply(lambda x: x.get('name'))
bridge = rd[rd['bridge'] == 'yes'].copy()
rd.drop(rd.index[rd['bridge'] == 'yes'], inplace = True)
rd = rd[rd['geometry'].type != 'Polygon']#) & (rd['lanes'] != np.nan)]
rd.dropna(subset=['lanes'], inplace=True)
#--- dissolve on street name
rd = rd.dissolve(by='name', as_index=False, dropna=False)
rd['length'] = rd.geometry.length
rd.sort_values('length', inplace=True, ascending=False)#, na_position='first')
rd.reset_index(drop=True, inplace=True)

# ---I want a vertex every 5 metre
from shapely.wkt import loads

def segmentize(geom):
    wkt = geom.wkt  # shapely Polygon to wkt
    geom = ogr.CreateGeometryFromWkt(wkt)  # create ogr geometry
    geom.Segmentize(5)  # densify geometry
    wkt2 = geom.ExportToWkt()  # ogr geometry to wkt
    new = loads(wkt2)  # wkt to shapely Polygon
    return new

rd['geometry'] = rd['geometry'].apply(segmentize)

#- buffer_func
def buffer(row):
    return row.geometry.buffer(round(float(row['lanes']) * 1.5, 2), 
                               cap_style=2, join_style=2)
rd_copy = rd.copy()
rd_copy['geometry'] = rd_copy.apply(buffer, axis=1)
rd_copy = rd_copy[rd_copy.geometry.length > 10]
rd_copy = rd_copy.explode()
rd_copy.reset_index(drop=True, inplace=True)
one = gpd.GeoDataFrame(rd_copy.iloc[0], crs='EPSG:32734').T

#--- no overlap
for i, row in rd_copy.iterrows():
    if i > 0:
        g = gpd.GeoDataFrame(row, crs='EPSG:32734').T
        g = gpd.overlay(g, one, how = 'difference')
        one = one.append(g)

#--- clean
one.reset_index(drop=True, inplace=True)
one = one.explode()
one.reset_index(drop=True, inplace=True

When I plot this:

fig, ax = plt.subplots(figsize=(15, 15))
ax.margins(-0.2, -0.1)
one.plot(ax=ax, facecolor='none', edgecolor='blue')

enter image description here

Notice; the traffic circle has artifacts. Why?
How do I clean this?
data here

2
  • 1
    Thanks for supplying code and data! There's a few things wrong with the code that errored when I cut and paste, but its close. Theres d_copy instead of rd_copy, and a bare EPSG:32734 which errors, and jparams from nowhere. Add the missing import lines and its perfect.
    – Spacedman
    Jun 7 at 16:39
  • fixed as per suggestions @Spacedman
    – arkriger
    Jun 8 at 9:40

1 Answer 1

4

Zooming in reveals...

enter image description here

...that the artefacts are due to butt-ends on line segments, because cap_style=2. When you do the buffer it creates flat ends at the end of each feature, so where you have two separate features you get these little gaps on the outside when you merge. The traffic circle is made of five features.

If you weren't trying to do per-feature buffer widths by number of lanes you would merge all these features into one and you'd get a smooth buffer, but that doesn't seem to be an option (although maybe there's a way of doing something by computing all features at all possible buffer widths and doing some operations but I've not thought this through). Or maybe merge all features with the same desired output buffer width and then buffer...

Possible cap_style=3 (square) will fix it - this makes the end-caps go past the end points by the width of the buffer, but you'll probably get little extra triangles where they overlap.

1
  • thank you @Spacedman
    – arkriger
    Jun 8 at 9:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.