CODEWITHSUNDEEP
pythonshapelysimplify
Paul Jurczak
Paul Jurczak

Reputation: 8173

How to get original indicies of a polygon after applying shapely.simplify?

This Python script:

from shapely import simplify, points, contains, Point

circle = Point(0, 0).buffer(1.0, quad_segs=8).exterior
simple = simplify(circle, 0.1)

simplifies polygon circle (red) and produces polygon simple (blue) with a subset of circle vertices:

enter image description here

iCircle contains the original indices of simple vertices: [0, 4, 8, 12, 16, 20, 24, 28, 32]:

iCircle = []

for i, p in enumerate(points(circle.coords)):
  if contains(simple, p):
    iCircle.append(i)

How to compute it without a costly lookup like the one above?

Circle is just a simple example. My question concerns an arbitrary polygon.

Upvotes: 0

Views: 43

Answers (2)

Pieter
Pieter

Reputation: 1516

It only supports numpy array lines as input, so you would need to do some plumbing to simplify shapely polygons, but you could check out the simplification library.

It's API supports returning the indices of the coordinates to retain after simplification instead of the end result... which sounds like what you are looking for.

I've used this library to make a simplify function where you can "blacklist" points on polygons from being simplified away, so for the plumbing part you could get some inspiration in pygeoops.simplify. A pull request to expose this feature in pygeoops is also an option...

Code sample:

import numpy as np
from simplification.cutil import simplify_coords_idx

# assume an array of coordinates: orig
simplified = simplify_coords_idx(orig, 1.0)
# build new geometry using only retained coordinates
orig_simplified = orig[simplified]

Upvotes: 2

Grismar
Grismar

Reputation: 31379

Given that Shapely uses the Douglas-Peucker algorithm for geometry smoothing, you're essentially asking how to determine what vertices it will retain for a certain distance parameter (epsilon), or 'tolerance' as Shapely calls it.

If you're only looking at a circle, every next vertex will at a constant index distance from the previous, due to the circle being perfectly regular. There's a decent explanation of the algorithm here, but for a circle, things are greatly simplified.

So, you could just find the indices for the first two, and the rest will follow.

But given that these will rarely be very large numbers, and Shapely uses numpy ndarray as a datatype for points(), you could just:

indices = np.where(np.in1d(points(circle.coords), points(simple.coords)))[0]
print(indices)

Output:

[ 0  4  8 12 16 20 24 28 32]

And that would work for other shapes you need this for as well.

Upvotes: 1

Related Questions