Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[Done] Fix WKB/WKT of empty (multi)points on GEOS 3.9 #392

Merged
merged 14 commits into from
Sep 18, 2021
Merged

[Done] Fix WKB/WKT of empty (multi)points on GEOS 3.9 #392

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
060738b
Some tests
caspervdw Sep 15, 2021
87f2f24
Fix
caspervdw Sep 15, 2021
6aeac42
Tests
caspervdw Sep 15, 2021
8edb2f6
Clean
caspervdw Sep 15, 2021
0fdbc1e
More tests and add GEOS 3.9.0 to matrix
caspervdw Sep 15, 2021
f1f7f19
Remove 3.9.0 test
caspervdw Sep 15, 2021
b6dd4b0
Remover conditional skip
caspervdw Sep 17, 2021
df262c1
Docs and tests
caspervdw Sep 17, 2021
e24743e
Remove the multipoint WKT empty fix
caspervdw Sep 17, 2021
55d2503
Fix test
caspervdw Sep 17, 2021
88d49ba
Update test_io.py
caspervdw Sep 17, 2021
2c792b5
Finish tests
caspervdw Sep 18, 2021
6fd61e4
Make test more readable
caspervdw Sep 18, 2021
e316c6f
Changelog [ci skip]
caspervdw Sep 18, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions pygeos/tests/common.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
"POLYGON((0 0, 0 10, 10 10, 10 0, 0 0), (2 2, 2 4, 4 4, 4 2, 2 2))"
)
empty_point = pygeos.Geometry("POINT EMPTY")
empty_point_z = pygeos.Geometry("POINT Z EMPTY")
empty_line_string = pygeos.Geometry("LINESTRING EMPTY")
empty_polygon = pygeos.Geometry("POLYGON EMPTY")
empty = pygeos.Geometry("GEOMETRYCOLLECTION EMPTY")
Expand Down
72 changes: 30 additions & 42 deletions pygeos/tests/test_io.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@

import pygeos

from .common import all_types, empty_point, point, point_z
from .common import all_types, empty_point, empty_point_z, point, point_z

# fmt: off
POINT11_WKB = b"\x01\x01\x00\x00\x00" + struct.pack("<2d", 1.0, 1.0)
Expand Down Expand Up @@ -375,17 +375,24 @@ def test_to_wkb_srid():


@pytest.mark.skipif(
pygeos.geos_version >= (3, 8, 0), reason="Pre GEOS 3.8.0 has 3D empty points"
pygeos.geos_version < (3, 9, 0),
reason="Pre GEOS 3.9 empty point serialization is buggy",
Copy link
Member

Choose a reason for hiding this comment

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

Just to ensure I am correctly reading this: you removed testing for this for GEOS <= 3.8? (which I am fine with, given the version dependent behaviour, as long as we can guarantee correct and stable behaviour starting with GEOS 3.9)

Copy link
Member

Choose a reason for hiding this comment

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

(although there is now a whole block of code in geos.c between #if !GEOS_SINCE_3_9_0 / #endif that isn't really tested at all? Or is it still tested through other tests?)

Copy link
Member Author

Choose a reason for hiding this comment

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

You’re right, there should be some testing in GEOS<3.9, now the code isn’t touched anymore in the unittests.

Not sure what exactly to test. I don’t like “confirming faulty behaviour” which the previous tests did. Maybe just run the function and make sure it doesn’t segfault. I’ll remove the skip to see what exactly we have to do on earlier versions.

)
@pytest.mark.parametrize(
"geom,dims,expected",
[
(empty_point, 2, POINT_NAN_WKB),
(empty_point, 3, POINTZ_NAN_WKB),
(empty_point, 3, POINT_NAN_WKB),
Copy link
Member

Choose a reason for hiding this comment

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

This behaviour changed? (which seems a correct change though)

Copy link
Member Author

Choose a reason for hiding this comment

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

On GEOS < 3.8, to_wkb(empty_point, output_dimension=3) indeed gave the (faulty) 3D WKB string
On GEOS ==3.8, this works correctly, but you always get a 2D WKB (also for empty_point_z
On GEOS >= 3.9, everything works fine.

I split the test into 3 to make this more clear and updated the docs.

(empty_point_z, 2, POINT_NAN_WKB),
(empty_point_z, 3, POINTZ_NAN_WKB),
(pygeos.multipoints([empty_point]), 2, MULTIPOINT_NAN_WKB),
(pygeos.multipoints([empty_point]), 3, MULTIPOINTZ_NAN_WKB),
(pygeos.multipoints([empty_point]), 3, MULTIPOINT_NAN_WKB),
(pygeos.multipoints([empty_point_z]), 2, MULTIPOINT_NAN_WKB),
(pygeos.multipoints([empty_point_z]), 3, MULTIPOINTZ_NAN_WKB),
(pygeos.geometrycollections([empty_point]), 2, GEOMETRYCOLLECTION_NAN_WKB),
(pygeos.geometrycollections([empty_point]), 3, GEOMETRYCOLLECTIONZ_NAN_WKB),
(pygeos.geometrycollections([empty_point]), 3, GEOMETRYCOLLECTION_NAN_WKB),
(pygeos.geometrycollections([empty_point_z]), 2, GEOMETRYCOLLECTION_NAN_WKB),
(pygeos.geometrycollections([empty_point_z]), 3, GEOMETRYCOLLECTIONZ_NAN_WKB),
(
pygeos.geometrycollections([pygeos.multipoints([empty_point])]),
2,
Expand All @@ -394,42 +401,21 @@ def test_to_wkb_srid():
(
pygeos.geometrycollections([pygeos.multipoints([empty_point])]),
3,
NESTED_COLLECTIONZ_NAN_WKB,
NESTED_COLLECTION_NAN_WKB,
),
],
)
def test_to_wkb_point_empty_pre_geos38(geom, dims, expected):
actual = pygeos.to_wkb(geom, output_dimension=dims, byte_order=1)
# Use numpy.isnan; there are many byte representations for NaN
assert actual[: -dims * 8] == expected[: -dims * 8]
assert np.isnan(struct.unpack("<{}d".format(dims), actual[-dims * 8 :])).all()


@pytest.mark.skipif(
pygeos.geos_version < (3, 8, 0), reason="Post GEOS 3.8.0 has 2D empty points"
)
@pytest.mark.parametrize(
"geom,dims,expected",
[
(empty_point, 2, POINT_NAN_WKB),
(empty_point, 3, POINT_NAN_WKB),
(pygeos.multipoints([empty_point]), 2, MULTIPOINT_NAN_WKB),
(pygeos.multipoints([empty_point]), 3, MULTIPOINT_NAN_WKB),
(pygeos.geometrycollections([empty_point]), 2, GEOMETRYCOLLECTION_NAN_WKB),
(pygeos.geometrycollections([empty_point]), 3, GEOMETRYCOLLECTION_NAN_WKB),
(
pygeos.geometrycollections([pygeos.multipoints([empty_point])]),
pygeos.geometrycollections([pygeos.multipoints([empty_point_z])]),
2,
NESTED_COLLECTION_NAN_WKB,
),
(
pygeos.geometrycollections([pygeos.multipoints([empty_point])]),
pygeos.geometrycollections([pygeos.multipoints([empty_point_z])]),
3,
NESTED_COLLECTION_NAN_WKB,
NESTED_COLLECTIONZ_NAN_WKB,
),
],
)
def test_to_wkb_point_empty_post_geos38(geom, dims, expected):
def test_to_wkb_point_empty(geom, dims, expected):
# Post GEOS 3.8: empty point is 2D
caspervdw marked this conversation as resolved.
Show resolved Hide resolved
actual = pygeos.to_wkb(geom, output_dimension=dims, byte_order=1)
# Use numpy.isnan; there are many byte representations for NaN
caspervdw marked this conversation as resolved.
Show resolved Hide resolved
Expand All @@ -438,24 +424,26 @@ def test_to_wkb_point_empty_post_geos38(geom, dims, expected):


@pytest.mark.parametrize(
"wkb,expected_type",
"wkb,expected_type,expected_dim",
[
(POINT_NAN_WKB, 0),
(POINTZ_NAN_WKB, 0),
(MULTIPOINT_NAN_WKB, 4),
(MULTIPOINTZ_NAN_WKB, 4),
(GEOMETRYCOLLECTION_NAN_WKB, 7),
(GEOMETRYCOLLECTIONZ_NAN_WKB, 7),
(NESTED_COLLECTION_NAN_WKB, 7),
(NESTED_COLLECTIONZ_NAN_WKB, 7),
(POINT_NAN_WKB, 0, 2),
(POINTZ_NAN_WKB, 0, 3),
(MULTIPOINT_NAN_WKB, 4, 2),
(MULTIPOINTZ_NAN_WKB, 4, 3),
(GEOMETRYCOLLECTION_NAN_WKB, 7, 2),
(GEOMETRYCOLLECTIONZ_NAN_WKB, 7, 3),
(NESTED_COLLECTION_NAN_WKB, 7, 2),
(NESTED_COLLECTIONZ_NAN_WKB, 7, 3),
],
)
def test_from_wkb_point_empty(wkb, expected_type):
def test_from_wkb_point_empty(wkb, expected_type, expected_dim):
geom = pygeos.from_wkb(wkb)
# POINT (nan nan) transforms to an empty point
# Note that the dimensionality (2D/3D) is GEOS-version dependent
assert pygeos.is_empty(geom)
assert pygeos.get_type_id(geom) == expected_type
# Note that the dimensionality (2D/3D) is GEOS-version dependent
if pygeos.geos_version >= (3, 9, 0):
assert pygeos.get_coordinate_dimension(geom) == expected_dim


def test_to_wkb_point_empty_srid():
Expand Down
6 changes: 3 additions & 3 deletions src/geos.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,8 @@ char multipoint_has_point_empty(GEOSContextHandle_t ctx, GEOSGeometry* geom) {
}

// POINT EMPTY is converted to POINT (nan nan)
// by GEOS >= 3.10.0. Before that, we do it ourselves here.
#if !GEOS_SINCE_3_10_0
// by GEOS >= 3.9.0. Before that, we do it ourselves here.
#if !GEOS_SINCE_3_9_0

/* Returns 1 if geometry is an empty point, 0 otherwise, 2 on error.
*/
Expand Down Expand Up @@ -255,7 +255,7 @@ GEOSGeometry* point_empty_to_nan_all_geoms(GEOSContextHandle_t ctx, GEOSGeometry
return result;
}

#endif // !GEOS_SINCE_3_10_0
#endif // !GEOS_SINCE_3_9_0

/* Checks whether the geometry is a multipoint with an empty point in it
*
Expand Down
4 changes: 2 additions & 2 deletions src/geos.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,11 +141,11 @@ extern PyObject* geos_exception[1];
extern void geos_error_handler(const char* message, void* userdata);
extern void geos_notice_handler(const char* message, void* userdata);
extern void destroy_geom_arr(void* context, GEOSGeometry** array, int length);
#if !GEOS_SINCE_3_10_0
#if !GEOS_SINCE_3_9_0
extern char has_point_empty(GEOSContextHandle_t ctx, GEOSGeometry* geom);
extern GEOSGeometry* point_empty_to_nan_all_geoms(GEOSContextHandle_t ctx,
GEOSGeometry* geom);
#endif // !GEOS_SINCE_3_10_0
#endif // !GEOS_SINCE_3_9_0
extern char check_to_wkt_compatible(GEOSContextHandle_t ctx, GEOSGeometry* geom);
extern char geos_interpolate_checker(GEOSContextHandle_t ctx, GEOSGeometry* geom);

Expand Down
6 changes: 3 additions & 3 deletions src/pygeom.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,8 +126,8 @@ static PyObject* GeometryObject_ToWKB(GeometryObject* obj) {

GEOS_INIT;

#if !GEOS_SINCE_3_10_0
// WKB Does not allow empty points in GEOS < 3.10.
#if !GEOS_SINCE_3_9_0
// WKB Does not allow empty points in GEOS < 3.9.
// We check for that and patch the POINT EMPTY if necessary
has_empty = has_point_empty(ctx, obj->ptr);
if (has_empty == 2) {
Expand All @@ -141,7 +141,7 @@ static PyObject* GeometryObject_ToWKB(GeometryObject* obj) {
}
#else
geom = obj->ptr;
#endif // !GEOS_SINCE_3_10_0
#endif // !GEOS_SINCE_3_9_0

/* Create the WKB writer */
writer = GEOSWKBWriter_create_r(ctx);
Expand Down
14 changes: 7 additions & 7 deletions src/ufuncs.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2791,9 +2791,9 @@ static void to_wkb_func(char** args, npy_intp* dimensions, npy_intp* steps, void
GEOSWKBWriter* writer;
unsigned char* wkb;
size_t size;
#if !GEOS_SINCE_3_10_0
#if !GEOS_SINCE_3_9_0
char has_empty;
#endif // !GEOS_SINCE_3_10_0
#endif // !GEOS_SINCE_3_9_0

if ((is2 != 0) || (is3 != 0) || (is4 != 0) || (is5 != 0)) {
PyErr_Format(PyExc_ValueError, "to_wkb function called with non-scalar parameters");
Expand Down Expand Up @@ -2835,8 +2835,8 @@ static void to_wkb_func(char** args, npy_intp* dimensions, npy_intp* steps, void
Py_INCREF(Py_None);
*out = Py_None;
} else {
#if !GEOS_SINCE_3_10_0
// WKB Does not allow empty points in GEOS<3.10.
#if !GEOS_SINCE_3_9_0
// WKB Does not allow empty points in GEOS<3.9.
// We check for that and patch the POINT EMPTY if necessary
has_empty = has_point_empty(ctx, in1);
if (has_empty == 2) {
Expand All @@ -2850,18 +2850,18 @@ static void to_wkb_func(char** args, npy_intp* dimensions, npy_intp* steps, void
}
#else
temp_geom = in1;
#endif // !GEOS_SINCE_3_10_0
#endif // !GEOS_SINCE_3_9_0
if (hex) {
wkb = GEOSWKBWriter_writeHEX_r(ctx, writer, temp_geom, &size);
} else {
wkb = GEOSWKBWriter_write_r(ctx, writer, temp_geom, &size);
}
#if !GEOS_SINCE_3_10_0
#if !GEOS_SINCE_3_9_0
// Destroy the temp_geom if it was patched (POINT EMPTY patch)
if (has_empty) {
GEOSGeom_destroy_r(ctx, temp_geom);
}
#endif // !GEOS_SINCE_3_10_0
#endif // !GEOS_SINCE_3_9_0
if (wkb == NULL) {
errstate = PGERR_GEOS_EXCEPTION;
goto finish;
Expand Down