5

I have this raster:

original raster

which I'd like to convert to vector data (linestring) (e.g. a Shapefile).

To this end, I'm using this Python code:

import rasterio
import geopandas as gpd
from skimage.morphology import skeletonize
from skimage.measure import label, regionprops

segmented_raster_path = '/path/to/segmented_raster.tif'

with rasterio.open(segmented_raster_path, "r") as src:
    img = src.read(1)
    transform = src.transform
    crs = src.crs

binary = (img > 0.1).astype(np.uint8) * 255
skel = skeletonize(binary).astype(np.uint8) * 255
labeled = label(skel, connectivity=2)
regions = regionprops(labeled)
lines = []

for region in regions:
    coords = region.coords  # row, col
    if len(coords) < 2:
        continue
    coords_xy = [transform * (col, row) for row, col in coords]  # convert to (x, y)
    if len(coords_xy) > 1:
        line = LineString(coords_xy)
        if line.length >= 1:  # filter very small elements
            lines.append(line)

gdf = gpd.GeoDataFrame(geometry=lines, crs=crs)
gdf.to_file('extracted_lines.shp')

This results in:

resulting raster

And if I zoom in one of the "large blue regions" in the middle-top, the line forms a zigzag pattern:

zooming in resulting raster

Why is the result so strange, and how can I get correct, smooth lines following the center of the white stripes?

1
  • 1
    Probably because when you find connected regions the resulting coordinates are defined in a row-major order. The algorithm has no concept of linearity Commented Oct 14 at 19:15

1 Answer 1

3

I don't know why this is happening.

An alternative could be to use pygeoops.centerline. Disclaimer: I'm the developer of pygeoops.

EDIT: Based on your comment, I think you would like that pixels that just touch at a corner are also treated as connected for skeleton/centerline generation. Hence, I added a buffer + dissolve + explode of the polygons before creating the centerlines which should accomplish this.

This gives this result (with default parameters):

Screenshot1

Zoomed in on the triangles you mentioned in your post, WITHOUT buffer + dissolve + explode, so pixels touching with a corner are not considered as connected:

Screenshot zoomed in, without buffer + dissolve

Zoomed in on the triangles you mentioned in your post, WITH buffer + dissolve + explode, so pixels just touching with a corner are considered as connected:

Screenshot zoomed in, with buffer + dissolve

Sample script including the buffer + dissolve + explode:

import geopandas as gpd
import matplotlib.pyplot as plt
import pygeoops
import shapely
import rasterio as rio
from rasterio import features

# Polygonize the image.
path = "https://i.sstatic.net/M5kSIlpB.png"
with rio.open(path) as src:
    image = src.read(1)

    # Flip image vertically to get the result like in the example.
    image = image[::-1, :]
    polygons = []
    for coords, value in features.shapes(image, image > 0):
        polygons.append(shapely.geometry.shape(coords))

    gdf = gpd.GeoDataFrame(geometry=polygons, crs=src.crs)

# Pixels/polygons that just touch at the corners should be considered connected as well
# for the centerline calculation. Hence, we buffer and dissolve the polygons first.
# Use a "mitre" join style for the buffer to avoid introducing too many points.
gdf["geometry"] = gdf.buffer(0.2, join_style="mitre")
gdf = gdf.dissolve().explode()

# Calculate the centerline.
centerlines_gdf = gdf.copy()
centerlines_gdf["geometry"] = pygeoops.centerline(gdf["geometry"])

print(centerlines_gdf)
centerlines_gdf.plot()
plt.show()
2
  • Thank your for your proposal, unfortunately with the geotiff the centerlines are split in many individual parts: i.sstatic.net/pzPwpitf.png Commented Oct 16 at 6:40
  • Based on the screenshot in your comment, I think you would like that pixels that only touch at a corner would also be treated as being connected for the skeleton/centerline generation... Hence, I added a buffer + dissolve + explode before creating the centerline which should accomplish this. Is that better? Commented Oct 16 at 14:05

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.