Nearest point functions now supports POINT Z

[jlucidar]

No description provided.

[coveralls]

Coverage increased (+0.02%) to 83.396% when pulling 46f5c84 on jlucidar:patch-1 into eae4dcb on Toblerity:master.

[sgillies]

@jlucidar can you point me to something in the GEOS change log or a GEOS PR or ticket discussing this?

[mwtoews]

Testing with GEOS master, it seems that GEOSNearestPoints_r only returns 2D geometries, regardless if the inputs are 3D. This change would need to happen in GEOS before any magic happens here.

[mwtoews]

I was about to close this, but checking again using:

from shapely import geos
from shapely.geometry import Point, LineString
from shapely.ops import nearest_points

print(geos.geos_version_string)

pt = Point(1, 2, 3)
ls = LineString([(0, 0, 20), (5, 0, 30)])
print([g.wkt for g in nearest_points(ls, pt)])
print([g.wkt for g in nearest_points(pt, ls)])

there appears to be some changes, e.g. with GEOS 3.9.1 these are strictly 2D:

3.9.1dev-CAPI-1.14.1
['POINT (1 0)', 'POINT (1 2)']
['POINT (1 2)', 'POINT (1 0)']

and with a development GEOS:

3.10.0dev-CAPI-1.15.0
['POINT (1 0)', 'POINT Z (1 2 3)']
['POINT Z (1 2 3)', 'POINT (1 0)']

I know GEOS is be capable of interpolating the Z coord for the other nearest point, because of this:

print(ls.interpolate(1.0).wkt)  # POINT Z (1 0 22)

xref https://trac.osgeo.org/geos/ticket/1102

[mwtoews]

But, also I've come to my senses that we should probably keep 2D geometry behaviour for Shapely 1.8 and earlier. Any future work on this end will not be via ctypes (as implemented in this PR).