Skip to content

Commit

Permalink
Change linearring closing logic (#431)
Browse files Browse the repository at this point in the history
  • Loading branch information
@caspervdw
caspervdw committed Nov 14, 2021
1 parent cfa4256 commit 4d1f33a
Show file tree
Hide file tree
Showing 6 changed files with 73 additions and 41 deletions.
4 changes: 3 additions & 1 deletion CHANGELOG.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,9 @@ Version 0.12 (unreleased)

**API Changes**

* ...
* When constructing a linearring through ``pygeos.linearrings`` or a polygon through
``pygeos.polygons`` the ring is automatically closed when supplied with 3 coordinates
also when the first and last are already equal (#431).

**Bug fixes**

Expand Down
11 changes: 7 additions & 4 deletions pygeos/_geometry.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,10 +137,13 @@ def simple_geometries_1d(object coordinates, object indices, int geometry_type,
# check if we need to close a linearring
if geometry_type == 2:
ring_closure = 0
for coord_idx in range(dims):
if coord_view[idx, coord_idx] != coord_view[idx + geom_size - 1, coord_idx]:
ring_closure = 1
break
if geom_size == 3:
ring_closure = 1
else:
for coord_idx in range(dims):
if coord_view[idx, coord_idx] != coord_view[idx + geom_size - 1, coord_idx]:
ring_closure = 1
break
# check the resulting size to prevent invalid rings
if geom_size + ring_closure < 4:
# the error equals PGERR_LINEARRING_NCOORDS (in pygeos/src/geos.h)
Expand Down
3 changes: 2 additions & 1 deletion pygeos/creation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,8 @@ def linestrings(coords, y=None, z=None, indices=None, out=None, **kwargs):
def linearrings(coords, y=None, z=None, indices=None, out=None, **kwargs):
"""Create an array of linearrings.
If the provided coords do not constitute a closed linestring, the first
If the provided coords do not constitute a closed linestring, or if there
are only 3 provided coords, the first
coordinate is duplicated at the end to close the ring. This function will
raise an exception if a linearring contains less than three points or if
the terminal coordinates contain NaN (not-a-number).
Expand Down
74 changes: 46 additions & 28 deletions pygeos/tests/test_creation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,20 +24,23 @@ def box_tpl(x1, y1, x2, y2):

def test_points_from_coords():
actual = pygeos.points([[0, 0], [2, 2]])
assert str(actual[0]) == "POINT (0 0)"
assert str(actual[1]) == "POINT (2 2)"
assert_geometries_equal(
actual, [pygeos.Geometry("POINT (0 0)"), pygeos.Geometry("POINT (2 2)")]
)


def test_points_from_xy():
actual = pygeos.points(2, [0, 1])
assert str(actual[0]) == "POINT (2 0)"
assert str(actual[1]) == "POINT (2 1)"
assert_geometries_equal(
actual, [pygeos.Geometry("POINT (2 0)"), pygeos.Geometry("POINT (2 1)")]
)


def test_points_from_xyz():
actual = pygeos.points(1, 1, [0, 1])
assert str(actual[0]) == "POINT Z (1 1 0)"
assert str(actual[1]) == "POINT Z (1 1 1)"
assert_geometries_equal(
actual, [pygeos.Geometry("POINT Z (1 1 0)"), pygeos.Geometry("POINT (1 1 1)")]
)


def test_points_invalid_ndim():
Expand All @@ -47,31 +50,42 @@ def test_points_invalid_ndim():

@pytest.mark.skipif(pygeos.geos_version < (3, 10, 0), reason="GEOS < 3.10")
def test_points_nan_becomes_empty():
assert str(pygeos.points(np.nan, np.nan)) == "POINT EMPTY"
actual = pygeos.points(np.nan, np.nan)
assert_geometries_equal(actual, pygeos.Geometry("POINT EMPTY"))


def test_linestrings_from_coords():
actual = pygeos.linestrings([[[0, 0], [1, 1]], [[0, 0], [2, 2]]])
assert str(actual[0]) == "LINESTRING (0 0, 1 1)"
assert str(actual[1]) == "LINESTRING (0 0, 2 2)"
assert_geometries_equal(
actual,
[
pygeos.Geometry("LINESTRING (0 0, 1 1)"),
pygeos.Geometry("LINESTRING (0 0, 2 2)"),
],
)


def test_linestrings_from_xy():
actual = pygeos.linestrings([0, 1], [2, 3])
assert str(actual) == "LINESTRING (0 2, 1 3)"
assert_geometries_equal(actual, pygeos.Geometry("LINESTRING (0 2, 1 3)"))


def test_linestrings_from_xy_broadcast():
x = [0, 1] # the same X coordinates for both linestrings
y = [2, 3], [4, 5] # each linestring has a different set of Y coordinates
actual = pygeos.linestrings(x, y)
assert str(actual[0]) == "LINESTRING (0 2, 1 3)"
assert str(actual[1]) == "LINESTRING (0 4, 1 5)"
assert_geometries_equal(
actual,
[
pygeos.Geometry("LINESTRING (0 2, 1 3)"),
pygeos.Geometry("LINESTRING (0 4, 1 5)"),
],
)


def test_linestrings_from_xyz():
actual = pygeos.linestrings([0, 1], [2, 3], 0)
assert str(actual) == "LINESTRING Z (0 2 0, 1 3 0)"
assert_geometries_equal(actual, pygeos.Geometry("LINESTRING Z (0 2 0, 1 3 0)"))


def test_linestrings_invalid_shape_scalar():
Expand All @@ -94,17 +108,26 @@ def test_linestrings_invalid_shape(shape):

def test_linearrings():
actual = pygeos.linearrings(box_tpl(0, 0, 1, 1))
assert str(actual) == "LINEARRING (1 0, 1 1, 0 1, 0 0, 1 0)"
assert_geometries_equal(
actual, pygeos.Geometry("LINEARRING (1 0, 1 1, 0 1, 0 0, 1 0)")
)


def test_linearrings_from_xy():
actual = pygeos.linearrings([0, 1, 2, 0], [3, 4, 5, 3])
assert str(actual) == "LINEARRING (0 3, 1 4, 2 5, 0 3)"
assert_geometries_equal(actual, pygeos.Geometry("LINEARRING (0 3, 1 4, 2 5, 0 3)"))


def test_linearrings_unclosed():
actual = pygeos.linearrings(box_tpl(0, 0, 1, 1)[:-1])
assert str(actual) == "LINEARRING (1 0, 1 1, 0 1, 0 0, 1 0)"
assert_geometries_equal(
actual, pygeos.Geometry("LINEARRING (1 0, 1 1, 0 1, 0 0, 1 0)")
)


def test_linearrings_unclosed_all_coords_equal():
actual = pygeos.linearrings([(0, 0), (0, 0), (0, 0)])
assert_geometries_equal(actual, pygeos.Geometry("LINEARRING (0 0, 0 0, 0 0, 0 0)"))


def test_linearrings_invalid_shape_scalar():
Expand All @@ -121,9 +144,6 @@ def test_linearrings_invalid_shape_scalar():
(2, 2, 2), # 2 linearrings of 2 2D points
(1, 2, 2), # 1 linearring of 2 2D points
(2, 2), # 1 linearring of 2 2D points (scalar)
(2, 3, 2), # 2 linearrings of 3 2D points
(1, 3, 2), # 1 linearring of 3 2D points
(3, 2), # 1 linearring of 3 2D points (scalar)
],
)
def test_linearrings_invalid_shape(shape):
Expand All @@ -145,7 +165,9 @@ def test_linearrings_all_nan():

def test_polygon_from_linearring():
actual = pygeos.polygons(pygeos.linearrings(box_tpl(0, 0, 1, 1)))
assert str(actual) == "POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))"
assert_geometries_equal(
actual, pygeos.Geometry("POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))")
)


def test_polygons_none():
Expand All @@ -155,7 +177,9 @@ def test_polygons_none():

def test_polygons():
actual = pygeos.polygons(box_tpl(0, 0, 1, 1))
assert str(actual) == "POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))"
assert_geometries_equal(
actual, pygeos.Geometry("POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))")
)


def test_polygon_no_hole_list_raises():
Expand Down Expand Up @@ -239,9 +263,6 @@ def test_polygons_not_enough_points_in_shell_scalar():
(2, 2, 2), # 2 linearrings of 2 2D points
(1, 2, 2), # 1 linearring of 2 2D points
(2, 2), # 1 linearring of 2 2D points (scalar)
(2, 3, 2), # 2 linearrings of 3 2D points
(1, 3, 2), # 1 linearring of 3 2D points
(3, 2), # 1 linearring of 3 2D points (scalar)
],
)
def test_polygons_not_enough_points_in_shell(shape):
Expand Down Expand Up @@ -269,9 +290,6 @@ def test_polygons_not_enough_points_in_holes_scalar():
(2, 2, 2), # 2 linearrings of 2 2D points
(1, 2, 2), # 1 linearring of 2 2D points
(2, 2), # 1 linearring of 2 2D points (scalar)
(2, 3, 2), # 2 linearrings of 3 2D points
(1, 3, 2), # 1 linearring of 3 2D points
(3, 2), # 1 linearring of 3 2D points (scalar)
],
)
def test_polygons_not_enough_points_in_holes(shape):
Expand All @@ -296,7 +314,7 @@ def test_polygons_not_enough_points_in_holes(shape):
)
def test_create_collection_only_none(func, expected):
actual = func(np.array([None], dtype=object))
assert str(actual) == expected
assert_geometries_equal(actual, pygeos.Geometry(expected))


@pytest.mark.parametrize(
Expand Down
6 changes: 5 additions & 1 deletion pygeos/tests/test_creation_indices.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,6 @@ def test_linearrings(coordinates):
@pytest.mark.parametrize(
"coordinates",
[
([[1, 1], [2, 1], [1, 1]]), # too few coordinates
([[1, np.nan], [2, 1], [2, 2], [1, 1]]), # starting with nan
],
)
Expand All @@ -198,6 +197,11 @@ def test_linearrings_invalid(coordinates):
pygeos.linearrings(coordinates, indices=np.zeros(len(coordinates)))


def test_linearrings_unclosed_all_coords_equal():
actual = pygeos.linearrings([(0, 0), (0, 0), (0, 0)], indices=np.zeros(3))
assert_geometries_equal(actual, pygeos.Geometry("LINEARRING (0 0, 0 0, 0 0, 0 0)"))


@pytest.mark.parametrize(
"indices,expected",
[
Expand Down
16 changes: 10 additions & 6 deletions src/ufuncs.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2252,12 +2252,16 @@ static void linearrings_func(char** args, npy_intp* dimensions, npy_intp* steps,
DOUBLE_COREDIM_LOOP_OUTER {
/* check if first and last coords are equal; duplicate if necessary */
ring_closure = 0;
DOUBLE_COREDIM_LOOP_INNER_2 {
first_coord = *(double*)(ip1 + i_c2 * cs2);
last_coord = *(double*)(ip1 + (n_c1 - 1) * cs1 + i_c2 * cs2);
if (first_coord != last_coord) {
ring_closure = 1;
break;
if (n_c1 == 3) {
ring_closure = 1;
} else {
DOUBLE_COREDIM_LOOP_INNER_2 {
first_coord = *(double*)(ip1 + i_c2 * cs2);
last_coord = *(double*)(ip1 + (n_c1 - 1) * cs1 + i_c2 * cs2);
if (first_coord != last_coord) {
ring_closure = 1;
break;
}
}
}
/* the minimum number of coordinates in a linearring is 4 */
Expand Down

0 comments on commit 4d1f33a

Please sign in to comment.