Skip to content

Commit

Permalink
Faster ring orientation checks, enforce geojson output ring orientation
Browse files Browse the repository at this point in the history
  • Loading branch information
@karimbahgat
karimbahgat committed Aug 31, 2021
1 parent 103bad0 commit 298d981
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 48 deletions.
65 changes: 49 additions & 16 deletions shapefile.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,14 +159,32 @@ class _Array(array.array):
def __repr__(self):
return str(self.tolist())

def signed_area(coords):
def signed_area(coords, fast=False):
"""Return the signed area enclosed by a ring using the linear time
algorithm. A value >= 0 indicates a counter-clockwise oriented ring.
A faster version is possible by setting 'fast' to True, which returns
2x the area, e.g. if you're only interested in the sign of the area.
"""
xs, ys = map(list, list(zip(*coords))[:2]) # ignore any z or m values
xs.append(xs[1])
ys.append(ys[1])
return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(coords)))/2.0
area2 = sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(coords)))
if fast:
return area2
else:
return area2 / 2.0

def is_cw(coords):
"""Returns True if a polygon ring has clockwise orientation, determined
by a negatively signed area.
"""
area2 = signed_area(coords, fast=True)
return area2 < 0

def rewind(coords):
"""Returns the input coords in reversed order.
"""
return list(reversed(coords))

def ring_bbox(coords):
"""Calculates and returns the bounding box of a ring.
Expand Down Expand Up @@ -260,7 +278,7 @@ def itercoords():
if not is_straight_line:
# get triplet orientation
closed_triplet = triplet + [triplet[0]]
triplet_ccw = signed_area(closed_triplet) >= 0
triplet_ccw = not is_cw(closed_triplet)
# check that triplet has the same orientation as the ring (means triangle is inside the ring)
if ccw == triplet_ccw:
# get triplet centroid
Expand Down Expand Up @@ -300,11 +318,13 @@ def organize_polygon_rings(rings, return_errors=None):
for ring in rings:
# shapefile format defines a polygon as a sequence of rings
# where exterior rings are clockwise, and holes counterclockwise
if signed_area(ring) < 0:
if is_cw(ring):
# ring is exterior
ring = rewind(ring) # GeoJSON and Shapefile exteriors have opposite orientation
exteriors.append(ring)
else:
# ring is a hole
ring = rewind(ring) # GeoJSON and Shapefile holes have opposite orientation
holes.append(ring)

# if only one exterior, then all holes belong to that exterior
Expand Down Expand Up @@ -340,7 +360,8 @@ def organize_polygon_rings(rings, return_errors=None):

if len(exterior_candidates) > 1:
# get hole sample point
hole_sample = ring_sample(holes[hole_i], ccw=True)
# Note: all rings now follow GeoJSON orientation, i.e. holes are clockwise
hole_sample = ring_sample(holes[hole_i], ccw=False)
# collect new exterior candidates
new_exterior_candidates = []
for ext_i in exterior_candidates:
Expand All @@ -357,7 +378,7 @@ def organize_polygon_rings(rings, return_errors=None):

if len(exterior_candidates) > 1:
# exterior candidate with the smallest area is the hole's most immediate parent
ext_i = sorted(exterior_candidates, key=lambda x: abs(signed_area(exteriors[x])))[0]
ext_i = sorted(exterior_candidates, key=lambda x: abs(signed_area(exteriors[x], fast=True)))[0]
hole_exteriors[hole_i] = [ext_i]

# separate out holes that are orphaned (not contained by any exterior)
Expand All @@ -383,7 +404,10 @@ def organize_polygon_rings(rings, return_errors=None):

# add orphan holes as exteriors
for hole_i in orphan_holes:
ext = holes[hole_i] # could potentially reverse their order, but in geojson winding order doesn't matter
ext = holes[hole_i]
# since this was previously a clockwise ordered hole, inverse the winding order
ext = rewind(ext)
# add as single exterior without any holes
poly = [ext]
polys.append(poly)

Expand All @@ -396,7 +420,10 @@ def organize_polygon_rings(rings, return_errors=None):
else:
if return_errors is not None:
return_errors['polygon_only_holes'] = len(holes)
exteriors = holes # could potentially reverse their order, but in geojson winding order doesn't matter
exteriors = holes
# since these were previously clockwise ordered holes, inverse the winding order
exteriors = [rewind(ext) for ext in exteriors]
# add as single exterior without any holes
polys = [[ext] for ext in exteriors]
return polys

Expand Down Expand Up @@ -577,12 +604,15 @@ def _from_geojson(geoj):
parts = []
index = 0
for i,ext_or_hole in enumerate(geoj["coordinates"]):
if i == 0 and not signed_area(ext_or_hole) < 0:
# although the latest GeoJSON spec states that exterior rings should have
# counter-clockwise orientation, we explicitly check orientation since older
# GeoJSONs might not enforce this.
if i == 0 and not is_cw(ext_or_hole):
# flip exterior direction
ext_or_hole = list(reversed(ext_or_hole))
elif i > 0 and not signed_area(ext_or_hole) >= 0:
ext_or_hole = rewind(ext_or_hole)
elif i > 0 and is_cw(ext_or_hole):
# flip hole direction
ext_or_hole = list(reversed(ext_or_hole))
ext_or_hole = rewind(ext_or_hole)
points.extend(ext_or_hole)
parts.append(index)
index += len(ext_or_hole)
Expand All @@ -604,12 +634,15 @@ def _from_geojson(geoj):
index = 0
for polygon in geoj["coordinates"]:
for i,ext_or_hole in enumerate(polygon):
if i == 0 and not signed_area(ext_or_hole) < 0:
# although the latest GeoJSON spec states that exterior rings should have
# counter-clockwise orientation, we explicitly check orientation since older
# GeoJSONs might not enforce this.
if i == 0 and not is_cw(ext_or_hole):
# flip exterior direction
ext_or_hole = list(reversed(ext_or_hole))
elif i > 0 and not signed_area(ext_or_hole) >= 0:
ext_or_hole = rewind(ext_or_hole)
elif i > 0 and is_cw(ext_or_hole):
# flip hole direction
ext_or_hole = list(reversed(ext_or_hole))
ext_or_hole = rewind(ext_or_hole)
points.extend(ext_or_hole)
parts.append(index)
index += len(ext_or_hole)
Expand Down
67 changes: 35 additions & 32 deletions test_shapefile.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@
],
[0],
{'type':'Polygon','coordinates':[
[(1,1),(1,9),(9,9),(9,1),(1,1)],
shapefile.rewind([(1,1),(1,9),(9,9),(9,1),(1,1)]),
]}
),
(shapefile.POLYGON, # single polygon, holes (ordered)
Expand All @@ -52,9 +52,9 @@
],
[0,5,5+5],
{'type':'Polygon','coordinates':[
[(1,1),(1,9),(9,9),(9,1),(1,1)], # exterior
[(2,2),(4,2),(4,4),(2,4),(2,2)], # hole 1
[(5,5),(7,5),(7,7),(5,7),(5,5)], # hole 2
shapefile.rewind([(1,1),(1,9),(9,9),(9,1),(1,1)]), # exterior
shapefile.rewind([(2,2),(4,2),(4,4),(2,4),(2,2)]), # hole 1
shapefile.rewind([(5,5),(7,5),(7,7),(5,7),(5,5)]), # hole 2
]}
),
(shapefile.POLYGON, # single polygon, holes (unordered)
Expand All @@ -65,9 +65,9 @@
],
[0,5,5+5],
{'type':'Polygon','coordinates':[
[(1,1),(1,9),(9,9),(9,1),(1,1)], # exterior
[(2,2),(4,2),(4,4),(2,4),(2,2)], # hole 1
[(5,5),(7,5),(7,7),(5,7),(5,5)], # hole 2
shapefile.rewind([(1,1),(1,9),(9,9),(9,1),(1,1)]), # exterior
shapefile.rewind([(2,2),(4,2),(4,4),(2,4),(2,2)]), # hole 1
shapefile.rewind([(5,5),(7,5),(7,7),(5,7),(5,5)]), # hole 2
]}
),
(shapefile.POLYGON, # multi polygon, no holes
Expand All @@ -77,10 +77,10 @@
[0,5],
{'type':'MultiPolygon','coordinates':[
[ # poly 1
[(1,1),(1,9),(9,9),(9,1),(1,1)],
shapefile.rewind([(1,1),(1,9),(9,9),(9,1),(1,1)]),
],
[ # poly 2
[(11,11),(11,19),(19,19),(19,11),(11,11)],
shapefile.rewind([(11,11),(11,19),(19,19),(19,11),(11,11)]),
],
]}
),
Expand All @@ -95,14 +95,14 @@
[0,5,10,15,20,25],
{'type':'MultiPolygon','coordinates':[
[ # poly 1
[(1,1),(1,9),(9,9),(9,1),(1,1)], # exterior
[(2,2),(4,2),(4,4),(2,4),(2,2)], # hole 1
[(5,5),(7,5),(7,7),(5,7),(5,5)], # hole 2
shapefile.rewind([(1,1),(1,9),(9,9),(9,1),(1,1)]), # exterior
shapefile.rewind([(2,2),(4,2),(4,4),(2,4),(2,2)]), # hole 1
shapefile.rewind([(5,5),(7,5),(7,7),(5,7),(5,5)]), # hole 2
],
[ # poly 2
[(11,11),(11,19),(19,19),(19,11),(11,11)], # exterior
[(12,12),(14,12),(14,14),(12,14),(12,12)], # hole 1
[(15,15),(17,15),(17,17),(15,17),(15,15)], # hole 2
shapefile.rewind([(11,11),(11,19),(19,19),(19,11),(11,11)]), # exterior
shapefile.rewind([(12,12),(14,12),(14,14),(12,14),(12,12)]), # hole 1
shapefile.rewind([(15,15),(17,15),(17,17),(15,17),(15,15)]), # hole 2
],
]}
),
Expand All @@ -116,15 +116,15 @@
[0,5,10,15,20],
{'type':'MultiPolygon','coordinates':[
[ # poly 1
[(1,1),(1,9),(9,9),(9,1),(1,1)], # exterior 1
[(2,2),(8,2),(8,8),(2,8),(2,2)], # hole 1.1
shapefile.rewind([(1,1),(1,9),(9,9),(9,1),(1,1)]), # exterior 1
shapefile.rewind([(2,2),(8,2),(8,8),(2,8),(2,2)]), # hole 1.1
],
[ # poly 2
[(3,3),(3,7),(7,7),(7,3),(3,3)], # exterior 2
[(4,4),(6,4),(6,6),(4,6),(4,4)], # hole 2.1
shapefile.rewind([(3,3),(3,7),(7,7),(7,3),(3,3)]), # exterior 2
shapefile.rewind([(4,4),(6,4),(6,6),(4,6),(4,4)]), # hole 2.1
],
[ # poly 3
[(4.5,4.5),(4.5,5.5),(5.5,5.5),(5.5,4.5),(4.5,4.5)], # exterior 3
shapefile.rewind([(4.5,4.5),(4.5,5.5),(5.5,5.5),(5.5,4.5),(4.5,4.5)]), # exterior 3
],
]}
),
Expand All @@ -138,15 +138,15 @@
[0,5,10,15,20+3],
{'type':'MultiPolygon','coordinates':[
[ # poly 1
[(1,1),(1,9),(9,9),(9,1),(1,1)], # exterior 1
[(2,2),(3,3),(4,2),(8,2),(8,8),(4,8),(2,8),(2,4),(2,2)], # hole 1.1
shapefile.rewind([(1,1),(1,9),(9,9),(9,1),(1,1)]), # exterior 1
shapefile.rewind([(2,2),(3,3),(4,2),(8,2),(8,8),(4,8),(2,8),(2,4),(2,2)]), # hole 1.1
],
[ # poly 2
[(3,3),(3,7),(7,7),(7,3),(3,3)], # exterior 2
[(4,4),(4,4),(6,4),(6,4),(6,4),(6,6),(4,6),(4,4)], # hole 2.1
shapefile.rewind([(3,3),(3,7),(7,7),(7,3),(3,3)]), # exterior 2
shapefile.rewind([(4,4),(4,4),(6,4),(6,4),(6,4),(6,6),(4,6),(4,4)]), # hole 2.1
],
[ # poly 3
[(4.5,4.5),(4.5,5.5),(5.5,5.5),(5.5,4.5),(4.5,4.5)], # exterior 3
shapefile.rewind([(4.5,4.5),(4.5,5.5),(5.5,5.5),(5.5,4.5),(4.5,4.5)]), # exterior 3
],
]}
),
Expand All @@ -162,16 +162,17 @@
[0,5,10,15,20,25,30],
{'type':'MultiPolygon','coordinates':[
[ # poly 1
[(1,1),(1,9),(9,9),(9,1),(1,1)], # exterior
[(2,2),(4,2),(4,4),(2,4),(2,2)], # hole 1
[(5,5),(7,5),(7,7),(5,7),(5,5)], # hole 2
shapefile.rewind([(1,1),(1,9),(9,9),(9,1),(1,1)]), # exterior
shapefile.rewind([(2,2),(4,2),(4,4),(2,4),(2,2)]), # hole 1
shapefile.rewind([(5,5),(7,5),(7,7),(5,7),(5,5)]), # hole 2
],
[ # poly 2
[(11,11),(11,19),(19,19),(19,11),(11,11)], # exterior
[(12,12),(14,12),(14,14),(12,14),(12,12)], # hole 1
[(15,15),(17,15),(17,17),(15,17),(15,15)], # hole 2
shapefile.rewind([(11,11),(11,19),(19,19),(19,11),(11,11)]), # exterior
shapefile.rewind([(12,12),(14,12),(14,14),(12,14),(12,12)]), # hole 1
shapefile.rewind([(15,15),(17,15),(17,17),(15,17),(15,15)]), # hole 2
],
[ # poly 3 (orphaned hole)
[ # poly 3 (orphaned hole)
# Note: due to the hole-to-exterior conversion, should return the same ring orientation
[(95,95),(97,95),(97,97),(95,97),(95,95)], # exterior
],
]}
Expand All @@ -183,9 +184,11 @@
[0,5],
{'type':'MultiPolygon','coordinates':[
[ # poly 1
# Note: due to the hole-to-exterior conversion, should return the same ring orientation
[(1,1),(9,1),(9,9),(1,9),(1,1)],
],
[ # poly 2
# Note: due to the hole-to-exterior conversion, should return the same ring orientation
[(11,11),(19,11),(19,19),(11,19),(11,11)],
],
]}
Expand Down

6 comments on commit 298d981

@QuLogic
Copy link

@QuLogic QuLogic commented on 298d981 Mar 6, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While the description says it's for enforcing GeoJSON output ring orientation, this seems to have an effect on reading shapefiles. Is that an intended change?

@karimbahgat
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right, that is actually what I meant, it enforces ring orientation in the output when creating a geojson, ie reading a shapefile to geojson. Previously I didn't do so, since geojson does not strictly require a particular ring orientation for outer/inner rings, it only recommends it. This commit enforced the recommended orientations.

@QuLogic
Copy link

@QuLogic QuLogic commented on 298d981 Mar 9, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I'm still a bit confused; I'm not using any GeoJSON. All I did was read a shapefile and look at it from the Python classes, and now the orientation is flipped. Are you saying pyshp classes are internally GeoJSON?

@karimbahgat
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The changes to ring orientation does not touch any of the code defining how a shapefile is read from file. Internally each shape is read into a Shape class with a type, a flat list of points and a flat list of parts. If you're reading the same shapefile before and after this commit and inspect the Shape contents they should be exactly the same and in the same order (eg sf.shape(0).points). The enforced ring orientation only applies when reading a shapefile and then converting from the internal shape structure to a geojson representation/dict (i.e. when calling the geo_interface property), or the opposite when writing a shapefile and adding a shape based on a geojson dict (i.e. sf.shape({'type':'Point', 'coordinates':[0,0]})). So I suspect the way you are "looking at" or inspecting the shape contents involves converting to geojson, eg via the __geo_interface__ attribute or some external library? The new code should at any rate be the "correct" orientation according to the geojson specification:

"Polygon rings MUST follow the right-hand rule for orientation (counterclockwise external rings, clockwise internal rings)".
https://datatracker.ietf.org/doc/html/rfc7946

Hope this makes sense. Are there any unforeseen consequences of the new enforcement of ring orientation for your application?

@QuLogic
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It breaks a test in Cartopy: SciTools/cartopy#2012 We don't directly use __geo_interface__, but I think this comes up implicitly via conversion to shapely?

It doesn't break any plots, so I think this is still safe. However, as noted in the issue, the result is now different from reading with Fiona, which could be annoying.

@karimbahgat
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that's likely what's happening, shapely using the geo interface to convert to shapely geometry. I added a detailed answer and possible solution in SciTools/cartopy#2012. In short, while this commit produces correct ring orientation according to the newest RFC7946 GeoJSON spec, the results are different because fiona/GDAL still use 2008 GeoJSON. After reading up on the RFC7946 specification in more detail, the required changes are so extensive that a complete shift to the new standard would be beyond the scope of pyshp. So the next version will likely switch back to not enforcing ring orientation in accordance with 2008 GeoJSON.

Please sign in to comment.