Reading/Creating a shapely.Point is 10s-100s times slower vs other Point-like objects

[idanmiara]

Expected behavior and actual behavior.

I've been doing some benchmarks (see below), which lead me to the conclusion that reading x/y/z from a shapely point is extremely slow, especially via the properties but also via get_coordinates.
Creating a shapely point is also extremely slow
(tens - hundreds times slower than the other containers)

Some rough speed factors compared to a plain tuple (may not be exactly from the output below as I multiple times):

read a single coordinate: tuple(1), Numpy(2.4), NamedTuple(1.14), DataClass(0.9) Shapely.Point(138)

reading 2 coords from a 2d point: tuple(1), Numpy(2.76), NamedTuple(1.14), DataClass(0.9), Shapely.Point(173 | get_coordinates: 45)

creating a 2d point: tuple(1), Numpy(15.4), NamedTuple(9.84), DataClass(7), ShapelyPoint(172-fastest method)

reading 3 coords from a 3d point: tuple(1), Numpy(3.43), NamedTuple(1.3), DataClass(0.9), Shapely.Point(364 | get_coordinates: 41)

creating a 3d point: tuple(1), Numpy(16.48), NamedTuple(10), DataClass(7), ShapelyPoint(173-fastest method)

Expected behavior and actual behavior.

Creating a Point and reading its coordinates should be reasonably fast.

Steps to reproduce the problem.

Run the following code:

from dataclasses import dataclass
from timeit import timeit
from typing import NamedTuple

import numpy as np
from shapely import Point as ShpPoint, points as shp_points
from shapely import get_coordinates


class TPoint2(NamedTuple):
    x: float
    y: float


class TPoint3(NamedTuple):
    x: float
    y: float
    z: float


@dataclass
class DPoint2:
    x: float
    y: float


@dataclass
class DPoint3:
    x: float
    y: float
    z: float


t2 = (1, 2)
tp2 = TPoint2(1, 2)
dp2 = DPoint2(1, 2)
sp2 = ShpPoint(1, 2)
np2 = np.array((1, 2))

t3 = (1, 2, 3)
tp3 = TPoint3(1, 2, 3)
dp3 = DPoint3(1, 2, 3)
sp3 = ShpPoint(1, 2, 3)
np3 = np.array((1, 2, 3))

iter = 1_000_000
for i in range(1):
    print(f'--- {i}: reading element from 2d Point')
    print(f'{timeit(lambda: t2[0], number=iter)=}')
    print(f'{timeit(lambda: np2[0], number=iter)=}')
    print(f'{timeit(lambda: tp2[0], number=iter)=}')
    print(f'{timeit(lambda: tp2.x, number=iter)=}')
    print(f'{timeit(lambda: dp2.x, number=iter)=}')
    print(f'{timeit(lambda: sp2.x, number=iter)=}')

    print(f'--- {i}: reading all from 2d Point')
    print(f'{timeit(lambda: (t2[0], t2[1]), number=iter)=}')
    print(f'{timeit(lambda: (np2[0], np2[1]), number=iter)=}')
    print(f'{timeit(lambda: tuple(np2), number=iter)=}')
    print(f'{timeit(lambda: (tp2[0], tp2[1]), number=iter)=}')
    print(f'{timeit(lambda: (tp2.x, tp2.y), number=iter)=}')
    print(f'{timeit(lambda: (dp2.x, dp2.y), number=iter)=}')
    print(f'{timeit(lambda: (sp2.x, sp2.y), number=iter)=}')
    print(f'{timeit(lambda: get_coordinates(sp2)[0], number=iter)=}')

    print(f'--- {i}: creating 2d Point')
    print(f'{timeit(lambda: (1, 2), number=iter)=}')
    print(f'{timeit(lambda: np.array((1, 2)), number=iter)=}')
    print(f'{timeit(lambda: TPoint2(1, 2), number=iter)=}')
    print(f'{timeit(lambda: DPoint2(1, 2), number=iter)=}')
    print(f'{timeit(lambda: ShpPoint(1, 2), number=iter)=}')
    print(f'{timeit(lambda: ShpPoint(1, 2), number=iter)=}')
    print(f'{timeit(lambda: shp_points(1, 2), number=iter)=}')
    print(f'{timeit(lambda: shp_points((1, 2)), number=iter)=}')

    print(f'--- {i}: reading element from 3d Point')
    print(f'{timeit(lambda: t3[2], number=iter)=}')
    print(f'{timeit(lambda: np3[2], number=iter)=}')
    print(f'{timeit(lambda: tp3[2], number=iter)=}')
    print(f'{timeit(lambda: tp3.z, number=iter)=}')
    print(f'{timeit(lambda: dp3.z, number=iter)=}')
    print(f'{timeit(lambda: sp3.z, number=iter)=}')

    print(f'--- {i}: reading all from 3d Point')
    print(f'{timeit(lambda: (t3[0], t3[1], t3[2]), number=iter)=}')
    print(f'{timeit(lambda: (np3[0], np3[1], np3[2]), number=iter)=}')
    print(f'{timeit(lambda: tuple(np3), number=iter)=}')
    print(f'{timeit(lambda: (tp3[0], tp3[1], tp3[2]), number=iter)=}')
    print(f'{timeit(lambda: (tp3.x, tp3.y, tp3.z), number=iter)=}')
    print(f'{timeit(lambda: (dp3.x, dp3.y, dp3.z), number=iter)=}')
    print(f'{timeit(lambda: (sp3.x, sp3.y, sp3.z), number=iter)=}')
    print(f'{timeit(lambda: get_coordinates(sp3, include_z=True)[0], number=iter)=}')

    print(f'--- {i}: creating 3d Point')
    print(f'{timeit(lambda: (1, 2, 3), number=iter)=}')
    print(f'{timeit(lambda: np.array((1, 2, 3)), number=iter)=}')
    print(f'{timeit(lambda: TPoint3(1, 2, 3), number=iter)=}')
    print(f'{timeit(lambda: DPoint3(1, 2, 3), number=iter)=}')
    print(f'{timeit(lambda: ShpPoint(1, 2, 3), number=iter)=}')
    print(f'{timeit(lambda: shp_points(1, 2, 3), number=iter)=}')
    print(f'{timeit(lambda: shp_points((1, 2, 3)), number=iter)=}')

My output was:

--- 0: reading element from 2d Point
timeit(lambda: t2[0], number=iter)=0.044886260002385825
timeit(lambda: np2[0], number=iter)=0.09956037401570939
timeit(lambda: tp2[0], number=iter)=0.05156597698805854
timeit(lambda: tp2.x, number=iter)=0.05421582600683905
timeit(lambda: dp2.x, number=iter)=0.03778560599312186
timeit(lambda: sp2.x, number=iter)=6.205369812989375
--- 0: reading all from 2d Point
timeit(lambda: (t2[0], t2[1]), number=iter)=0.06794292901759036
timeit(lambda: (np2[0], np2[1]), number=iter)=0.1881827589822933
timeit(lambda: tuple(np2), number=iter)=0.5759876949887257
timeit(lambda: (tp2[0], tp2[1]), number=iter)=0.08532903200830333
timeit(lambda: (tp2.x, tp2.y), number=iter)=0.0910857769777067
timeit(lambda: (dp2.x, dp2.y), number=iter)=0.06309307098854333
timeit(lambda: (sp2.x, sp2.y), number=iter)=11.755082368996227
timeit(lambda: get_coordinates(sp2)[0], number=iter)=3.1168896360031795
--- 0: creating 2d Point
timeit(lambda: (1, 2), number=iter)=0.03463764698244631
timeit(lambda: np.array((1, 2)), number=iter)=0.5365576750191394
timeit(lambda: TPoint2(1, 2), number=iter)=0.32849500098382123
timeit(lambda: DPoint2(1, 2), number=iter)=0.22573559000738896
timeit(lambda: ShpPoint(1, 2), number=iter)=6.568336394004291
timeit(lambda: ShpPoint(1, 2), number=iter)=6.4917032179946546
timeit(lambda: shp_points(1, 2), number=iter)=11.928747585014207
timeit(lambda: shp_points((1, 2)), number=iter)=5.85127081401879
--- 0: reading element from 3d Point
timeit(lambda: t3[2], number=iter)=0.04005866500665434
timeit(lambda: np3[2], number=iter)=0.09451862299465574
timeit(lambda: tp3[2], number=iter)=0.050292653002543375
timeit(lambda: tp3.z, number=iter)=0.05330886802403256
timeit(lambda: dp3.z, number=iter)=0.041282758000306785
timeit(lambda: sp3.z, number=iter)=16.58963368300465
--- 0: reading all from 3d Point
timeit(lambda: (t3[0], t3[1], t3[2]), number=iter)=0.07935839900164865
timeit(lambda: (np3[0], np3[1], np3[2]), number=iter)=0.2724746929889079
timeit(lambda: tuple(np3), number=iter)=0.616603036003653
timeit(lambda: (tp3[0], tp3[1], tp3[2]), number=iter)=0.10378881398355588
timeit(lambda: (tp3.x, tp3.y, tp3.z), number=iter)=0.11160819401266053
timeit(lambda: (dp3.x, dp3.y, dp3.z), number=iter)=0.07096943099168129
timeit(lambda: (sp3.x, sp3.y, sp3.z), number=iter)=28.769336715020472
timeit(lambda: get_coordinates(sp3, include_z=True)[0], number=iter)=3.2244306430220604
--- 0: creating 3d Point
timeit(lambda: (1, 2, 3), number=iter)=0.033816418988863006
timeit(lambda: np.array((1, 2, 3)), number=iter)=0.5440691379772034
timeit(lambda: TPoint3(1, 2, 3), number=iter)=0.33719124301569536
timeit(lambda: DPoint3(1, 2, 3), number=iter)=0.23762070699012838
timeit(lambda: ShpPoint(1, 2, 3), number=iter)=6.520570991007844
timeit(lambda: shp_points(1, 2, 3), number=iter)=13.25042435200885
timeit(lambda: shp_points((1, 2, 3)), number=iter)=5.8881938550039195

Operating system

Mac OS X 13.4

Shapely version and provenance

2.0.1 installed from PyPI using pip

[mwtoews]

A bit of déjà vu from #983 (comment) and following comments.

I'm unsure if a "single-point constructor" was ever written that optimizes for these cases. Users are encouraged to create an array of (e.g.) 1_000_000 point geometries, rather than repeating single point geometry operations.

[idanmiara]

A bit of déjà vu from #983 (comment) and following comments.

I'm unsure if a "single-point constructor" was ever written that optimizes for these cases. Users are encouraged to create an array of (e.g.) 1_000_000 point geometries, rather than repeating single point geometry operations.

Thanks for the reference.
Optimizing for arrays is good, but single object performance should also be reasonable IMOH.
Maybe an approach like pyproj4/pyproj#1206 can be applied here?

Also if we provide [x,y,z] properties, we somewhat encourage people to use them.
If the performance issue can be resolved at least for reading the coordinate values, maybe we can also add __get_item__ to allow to access p[0], p[1], p[2] as an alternative to p.x, p.y, p.z in the spirit of the same duality of NamedTuple point in the example above.

[jorisvandenbossche]

A bit of déjà vu from #983 (comment) and following comments.

Related to this slowdown in scalar operations compared to 1.8 (because of the refactor using array ufuncs for everything), and specifically for the Point constructor example discussed here: I did some effort to speed up the scalar Point(..) constructor(#1547), and currently this specific operation is actually slightly faster compared to 1.8:

In [5]: %timeit Point(0.0, 0.0)
7.36 µs ± 136 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)   # shapely 1.8.5.post1
5.21 µs ± 408 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)   # shapely 2.0.1

Now, that's 1) only one specific case, and 2) it's not because it does OK compared to 1.8 that it can't still be faster.

Optimizing for arrays is good, but single object performance should also be reasonable IMOH.
Maybe an approach like pyproj4/pyproj#1206 can be applied here?

I am not sure that is doable in general, because if I understand that PR correctly, it is essentially creating a special case implementation for the scalar case. But in shapely we have more than 100 ufuncs binding some GEOS method that all could have a separate scalar version.

Specifically for the Point(..) creation, this is certainly doable: we actually already have a scalar C version create_point that could be exposed as a Python function to use in the Point(..) constructor. And it might be worth doing that for the Point constructor, but in general, the question is where to stop.

(and I know you actually propose to stop at just providing coordinate access for Points, but what if then someone else passes by requesting this for one specific other operation?)

Also if we provide [x,y,z] properties, we somewhat encourage people to use them.

In my mind, it is encouraged to use those, if you are working with scalar objects. And if you have a numpy array of points, then the attribute also doesn't work out of the box (you have to write a for loop), and in that case the ufunc is both easier to use and faster.

I think in many applications the performance of those attributes won't be a problem (and if it is, it might be worth looking into the vectorized functions anyway, even if we can make those attributes a bit faster).

And to be clear, that's not to say that we shouldn't make the Geometry object attributes faster, if that is easy to do we for sure should do it. But it is a tradeoff with the added code complexity and maintenance to do so.

maybe we can also add __get_item__ to allow to access p[0], p[1], p[2] as an alternative to p.x, p.y, p.z in the spirit of the same duality of NamedTuple point in the example above.

In general, the fact that geometry objects don't have sequence-like behaviours is intentional: this was one of the design changes of shapely 2.0, and we actually provided this feature in 1.x (well, not exactly for Points, but we did for multi-geometries, and we did provide the array-interface for Points).
The idea now is that Point/LineString objects are not like tuples or sequences of coordinates, it's a single scalar object, that just happens to be backed by a set of coordinates. Treating geometries as sequence-like objects gave all kinds of corner cases when working with arrays of such objects. See https://shapely.readthedocs.io/en/stable/migration.html#multi-part-geometries-will-no-longer-be-sequences-length-iterable-indexable and the section afterwards about numpy array interface.

[idanmiara]

@jorisvandenbossche Thanks for your comments.
Good point about the unpacking issue if adding __get_item__.

[gatopeich]

Profiling the code with Py-spy, I see most CPU is spent on "wrapped" decorator:

This code is checking if a box Polygon contains each point of a long list of coordinates, like the last line in this snippet:

class Area:
    def __init__(self, xyz1, xyz2):
            self.polygon = shapely.polygon_box(*xyz1, *xyz2)

    def __contains__(self, x, y, z):
        return self.polygon.contains(Point(x, y, z))