Comparing cuspatial buffers vs Arrow buffers¶

For context see https://github.com/geopandas/geo-arrow-spec/issues/4

Test data¶

1
import geopandas

1
df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))

1
df
pop_est continent name iso_a3 gdp_md_est geometry
0 920938 Oceania Fiji FJI 8374.0 MULTIPOLYGON (((180.00000 -16.06713, 180.00000...
1 53950935 Africa Tanzania TZA 150600.0 POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...
2 603253 Africa W. Sahara ESH 906.5 POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3 35623680 North America Canada CAN 1674000.0 MULTIPOLYGON (((-122.84000 49.00000, -122.9742...
4 326625791 North America United States of America USA 18560000.0 MULTIPOLYGON (((-122.84000 49.00000, -120.0000...
... ... ... ... ... ... ...
172 7111024 Europe Serbia SRB 101800.0 POLYGON ((18.82982 45.90887, 18.82984 45.90888...
173 642550 Europe Montenegro MNE 10610.0 POLYGON ((20.07070 42.58863, 19.80161 42.50009...
174 1895250 Europe Kosovo -99 18490.0 POLYGON ((20.59025 41.85541, 20.52295 42.21787...
175 1218208 North America Trinidad and Tobago TTO 43570.0 POLYGON ((-61.68000 10.76000, -61.10500 10.890...
176 13026129 Africa S. Sudan SSD 20880.0 POLYGON ((30.83385 3.50917, 29.95350 4.17370, ...

177 rows × 6 columns

Conversion to Arrow nested ListArray / buffers¶

A quick conversion by extracting the GeoJSON-like nested list of coordinates for each of the geometries, and using the PyArrow pa.array() conversion utility to convert this into a ListArray, from which the underlying buffers can be extracted:

1
import pyarrow as pa

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
from shapely.geometry import MultiPolygon, Polygon

def get_nested_coordinates_multipolygon(data):
    
    nested_coordinates = []
    
    for geom in data:
        if isinstance(geom, MultiPolygon):
            nested_coordinates.append(geom.__geo_interface__["coordinates"])
        elif isinstance(geom, Polygon):
            nested_coordinates.append([geom.__geo_interface__["coordinates"]])
        else:
            raise TypeError(type(geom))

    return nested_coordinates

1
nested_coordinates = get_nested_coordinates_multipolygon(df.geometry)

1
2
3
4
5
6
multipolygon_type = pa.list_(
    pa.field("parts", pa.list_(
        pa.field("rings", pa.list_(
            pa.field("vertices", pa.list_(
                pa.field("xy", pa.float64(), nullable=False), 2),
            nullable=False)), nullable=False))))

1
arr = pa.array(nested_coordinates, multipolygon_type)

1
2
coords = np.asarray(arr.values.values.values.values)
coords

array([180.        , -16.06713266, 180.        , ...,   3.7819    ,
        30.83385242,   3.5091716 ])
1
offsets_multipolygons = np.asarray(arr.offsets)

1
offsets_polygons = np.asarray(arr.values.offsets)

1
offsets_rings = np.asarray(arr.values.values.offsets)

Conversion to cuspatial GeoArrow buffers¶

The code to run this is included at the bottom of the notebook as a copy from the cuspatial source code. Using this allows us to get the buffers cuspatial creates from GeoPandas data (which are then used to pass to the GPU memory):

1
adapter = GeoPandasAdapter(df.geometry)

1
2
meta = adapter.get_geopandas_meta()
meta.keys()

dict_keys(['input_types', 'input_lengths', 'inputs'])
1
2
buffers = adapter.get_geoarrow_host_buffers()
buffers.keys()

dict_keys(['points_xy', 'mpoints_xy', 'mpoints_offsets', 'lines_xy', 'lines_offsets', 'mlines', 'polygons_xy', 'polygons_polygons', 'polygons_rings', 'mpolygons'])

Comparing the results¶

The coordinates are exactly the same (interleaved xy array):

1
np.allclose(buffers["polygons_xy"], coords)

True

The offsets into the coordinates to indicate the start of rings is equivalent, only cuspatial uses actual offsets into the coordinate array, while for Arrow you need to multiply this with 2 (this is because for Arrow the inner field is a fixed-size list array (size 2), and the offsets are relative to the logical items (lists of 2) of this list array, not relative to the physical coordinates array. But because this is a fixed-size list array, we know the conversion is a simple * 2):

1
np.allclose(buffers["polygons_rings"], offsets_rings * 2)

True

The offsets where each Polygon starts (exterior ring, potential interior rings) in the ring offsets is exactly the same:

1
np.allclose(buffers["polygons_polygons"], offsets_polygons)

True

Finally, there are offsets / indices where each MultiPolygon starts (consisting of multiple Polygons).

Here, there are some differences because cuspatial keeps more closely track of which geometries are Polygon and which are MultiPolygon. While the Arrow list array simply assumes a len-1 MultiPolygon for the Polygons.

For Arrow this gives:

1
offsets_multipolygons

array([  0,   3,   4,   5,  35,  45,  46,  47,  51,  64,  66,  68,  69,
        70,  71,  72,  73,  74,  75,  89,  92,  93,  97,  98,  99, 100,
       101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113,
       114, 115, 116, 117, 118, 121, 122, 123, 124, 125, 126, 127, 128,
       129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
       142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 153, 154, 155,
       156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 168, 170,
       171, 172, 173, 174, 175, 177, 178, 179, 180, 181, 182, 183, 184,
       185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197,
       198, 199, 200, 201, 202, 203, 204, 206, 208, 209, 210, 211, 212,
       213, 214, 215, 216, 217, 218, 223, 225, 227, 228, 230, 231, 234,
       236, 238, 239, 241, 242, 249, 251, 252, 253, 254, 255, 256, 257,
       260, 261, 262, 263, 271, 272, 273, 274, 275, 276, 277, 278, 279,
       280, 281, 282, 283, 284, 285, 286, 287, 288], dtype=int32)

While cuspatial has only a subset of those values in the "mpolygons" buffer (I suppose this is only for the MultiPolygons in the array?):

1
print(buffers["mpolygons"])

[0, 3, 5, 35, 35, 45, 47, 51, 51, 64, 64, 66, 66, 68, 75, 89, 89, 92, 93, 97, 118, 121, 151, 153, 166, 168, 168, 170, 175, 177, 204, 206, 206, 208, 218, 223, 223, 225, 225, 227, 228, 230, 231, 234, 234, 236, 236, 238, 239, 241, 242, 249, 249, 251, 257, 260, 263, 271]

However, cuspatial also keeps track of the length of each input feature in the metadata:

1
print(meta["input_lengths"])

[3, 1, 1, 30, 10, 1, 1, 4, 13, 2, 2, 1, 1, 1, 1, 1, 1, 1, 14, 3, 1, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 5, 2, 2, 1, 2, 1, 3, 2, 2, 1, 2, 1, 7, 2, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

And based on this you can actually get the same offsets:

1
np.insert(np.cumsum(meta["input_lengths"]), 0, 0)

array([  0,   3,   4,   5,  35,  45,  46,  47,  51,  64,  66,  68,  69,
        70,  71,  72,  73,  74,  75,  89,  92,  93,  97,  98,  99, 100,
       101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113,
       114, 115, 116, 117, 118, 121, 122, 123, 124, 125, 126, 127, 128,
       129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
       142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 153, 154, 155,
       156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 168, 170,
       171, 172, 173, 174, 175, 177, 178, 179, 180, 181, 182, 183, 184,
       185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197,
       198, 199, 200, 201, 202, 203, 204, 206, 208, 209, 210, 211, 212,
       213, 214, 215, 216, 217, 218, 223, 225, 227, 228, 230, 231, 234,
       236, 238, 239, 241, 242, 249, 251, 252, 253, 254, 255, 256, 257,
       260, 261, 262, 263, 271, 272, 273, 274, 275, 276, 277, 278, 279,
       280, 281, 282, 283, 284, 285, 286, 287, 288])
1
np.allclose(np.insert(np.cumsum(meta["input_lengths"]), 0, 0), offsets_multipolygons)

True

The buffers["mpolygons"] are indeed the start / and offsets for only the MultiPolygons:

1
print(list(offsets_multipolygons[:-1][np.array(meta["input_lengths"]) > 1]))

[0, 5, 35, 47, 51, 64, 66, 75, 89, 93, 118, 151, 166, 168, 175, 204, 206, 218, 223, 225, 228, 231, 234, 236, 239, 242, 249, 257, 263]

1
print(buffers["mpolygons"][::2])

[0, 5, 35, 47, 51, 64, 66, 75, 89, 93, 118, 151, 166, 168, 175, 204, 206, 218, 223, 225, 228, 231, 234, 236, 239, 242, 249, 257, 263]

Conversion code from cuspatial¶

Copied from https://github.com/rapidsai/cuspatial/blob/main/python/cuspatial/cuspatial/io/geopandas_adapter.py

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
# Copyright (c) 2020-2021 NVIDIA CORPORATION.

import numpy as np
from geopandas import GeoSeries as gpGeoSeries
from shapely.geometry import (
    LineString,
    MultiLineString,
    MultiPoint,
    MultiPolygon,
    Point,
    Polygon,
)


class GeoPandasAdapter:
    def __init__(self, geoseries: gpGeoSeries):
        """
        GeoPandasAdapter copies a GeoPandas GeoSeries object iteratively into
        a set of arrays: points, multipoints, lines, and polygons.
        Parameters
        ----------
        geoseries : A GeoPandas GeoSeries
        """
        self.offsets = self._load_geometry_offsets(geoseries)
        self.buffers = self._read_geometries(geoseries, self.offsets)

    def _load_geometry_offsets(self, geoseries: gpGeoSeries) -> dict:
        """
        Computes the offet arrays and buffer sizes  that will be required
        to store the geometries.
        Parameters
        ----------
        geoseries : A GeoPandas GeoSeries
        """
        offsets = {
            "points": [0],
            "multipoints": [0],
            "lines": [0],
            "mlines": [],
            "polygons": {"polygons": [0], "rings": [0], "mpolys": []},
        }
        for geometry in geoseries:
            if isinstance(geometry, Point):
                # a single Point geometry will go into the GpuPoints
                # structure. No offsets are required, but an index to the
                # position in the GeoSeries is required.
                current = offsets["points"][-1]
                offsets["points"].append(len(geometry.xy) + current)
            elif isinstance(geometry, MultiPoint):
                # A MultiPoint geometry also is copied into the GpuPoints
                # structure. A MultiPoint object must be created, containing
                # the size of the number of points, the position they are
                # stored in GpuPoints, and the index of the MultiPoint in the
                # GeoSeries.
                current = offsets["multipoints"][-1]
                offsets["multipoints"].append(len(geometry) * 2 + current)
            elif isinstance(geometry, LineString):
                # A LineString geometry is stored in the GpuLines structure.
                # Every LineString has a size which is stored in the GpuLines
                # structure. The index of the LineString back into the
                # GeoSeries is also stored.
                current = offsets["lines"][-1]
                offsets["lines"].append(2 * len(geometry.coords) + current)
            elif isinstance(geometry, MultiLineString):
                # A MultiLineString geometry is stored identically to
                # LineString in the GpuLines structure. The index of the
                # GeoSeries object is also stored.
                offsets["mlines"].append(len(offsets["lines"]) - 1)
                for linestring in geometry:
                    current = offsets["lines"][-1]
                    offsets["lines"].append(
                        2 * len(linestring.coords) + current
                    )
                offsets["mlines"].append(len(offsets["lines"]) - 1)
            elif isinstance(geometry, Polygon):
                # A Polygon geometry is stored like a LineString and also
                # contains a buffer of sizes for each inner ring.
                num_rings = 1
                rings_current = offsets["polygons"]["rings"][-1]
                offsets["polygons"]["rings"].append(
                    len(geometry.exterior.coords) * 2 + rings_current
                )
                for interior in geometry.interiors:
                    rings_current = offsets["polygons"]["rings"][-1]
                    offsets["polygons"]["rings"].append(
                        len(interior.coords) * 2 + rings_current
                    )
                    num_rings = num_rings + 1
                current = offsets["polygons"]["polygons"][-1]
                offsets["polygons"]["polygons"].append(num_rings + current)
            elif isinstance(geometry, MultiPolygon):
                current = offsets["polygons"]["polygons"][-1]
                offsets["polygons"]["mpolys"].append(
                    len(offsets["polygons"]["polygons"]) - 1
                )
                for poly in geometry:
                    current = offsets["polygons"]["polygons"][-1]
                    num_rings = 1
                    rings_current = offsets["polygons"]["rings"][-1]
                    offsets["polygons"]["rings"].append(
                        len(poly.exterior.coords) * 2 + rings_current
                    )
                    for interior in poly.interiors:
                        rings_current = offsets["polygons"]["rings"][-1]
                        offsets["polygons"]["rings"].append(
                            len(interior.coords) * 2 + rings_current
                        )
                        num_rings = num_rings + 1
                    offsets["polygons"]["polygons"].append(num_rings + current)
                offsets["polygons"]["mpolys"].append(
                    len(offsets["polygons"]["polygons"]) - 1
                )
        return offsets

    def _read_geometries(self, geoseries: gpGeoSeries, offsets: dict,) -> dict:
        """
        Creates a set of buffers sized to fit all of the geometries and
        iteratively populates them with geometry coordinate values.
        Parameters
        ----------
        geoseries : A GeoPandas GeoSeries object.
        offsets : The set of offsets that correspond to the geoseries argument.
        """
        buffers = {
            "points": np.zeros(offsets["points"][-1]),
            "multipoints": np.zeros(offsets["multipoints"][-1]),
            "lines": np.zeros(offsets["lines"][-1]),
            "polygons": {
                "polygons": np.zeros(len(offsets["polygons"]["polygons"])),
                "rings": np.zeros(len(offsets["polygons"]["rings"])),
                "coords": np.zeros(offsets["polygons"]["rings"][-1]),
            },
        }
        read_count = {
            "points": 0,
            "multipoints": 0,
            "lines": 0,
            "polygons": 0,
        }
        inputs = []
        input_types = []
        input_lengths = []
        for geometry in geoseries:
            if isinstance(geometry, Point):
                # write a point to the points buffer
                # increase read_count of points pass
                i = read_count["points"] * 2
                buffers["points"][i] = geometry.x
                buffers["points"][i + 1] = geometry.y
                read_count["points"] = read_count["points"] + 1
                input_types.append("p")
                input_lengths.append(1)
                inputs.append({"type": "p", "length": 1})
            elif isinstance(geometry, MultiPoint):
                points = np.array(geometry)
                size = points.shape[0] * 2
                i = read_count["multipoints"]
                buffers["multipoints"][slice(i, i + size, 2)] = points[:, 0]
                buffers["multipoints"][slice(i + 1, i + size, 2)] = points[
                    :, 1
                ]
                read_count["multipoints"] = read_count["multipoints"] + size
                input_types.append("mp")
                input_lengths.append(len(geometry))
                inputs.append({"type": "mp", "length": len(geometry)})
            elif isinstance(geometry, LineString):
                size = len(geometry.xy[0]) * 2
                i = read_count["lines"]
                buffers["lines"][slice(i, i + size, 2)] = geometry.xy[0]
                buffers["lines"][slice(i + 1, i + size, 2)] = geometry.xy[1]
                read_count["lines"] = read_count["lines"] + size
                input_types.append("l")
                input_lengths.append(1)
                inputs.append({"type": "l", "length": 1})
            elif isinstance(geometry, MultiLineString):
                substrings = []
                for linestring in geometry:
                    size = len(linestring.xy[0]) * 2
                    i = read_count["lines"]
                    buffers["lines"][slice(i, i + size, 2)] = linestring.xy[0]
                    buffers["lines"][
                        slice(i + 1, i + size, 2)
                    ] = linestring.xy[1]
                    read_count["lines"] = read_count["lines"] + size
                    substrings.append({"type": "l", "length": size})
                input_types.append("ml")
                input_lengths.append(len(geometry))
                inputs.append(
                    {
                        "type": "ml",
                        "length": len(geometry),
                        "children": substrings,
                    }
                )
            elif isinstance(geometry, Polygon):
                # copy exterior
                exterior = geometry.exterior.coords.xy
                size = len(exterior[0]) * 2
                i = read_count["polygons"]
                buffers["polygons"]["coords"][
                    slice(i, i + size, 2)
                ] = exterior[0]
                buffers["polygons"]["coords"][
                    slice(i + 1, i + size, 2)
                ] = exterior[1]
                read_count["polygons"] = read_count["polygons"] + size
                interiors = geometry.interiors
                for interior in interiors:
                    interior_coords = interior.coords.xy
                    size = len(interior_coords[0]) * 2
                    i = read_count["polygons"]
                    buffers["polygons"]["coords"][
                        slice(i, i + size, 2)
                    ] = interior_coords[0]
                    buffers["polygons"]["coords"][
                        slice(i + 1, i + size, 2)
                    ] = interior_coords[1]
                    read_count["polygons"] = read_count["polygons"] + size
                input_types.append("poly")
                input_lengths.append(1)
                inputs.append({"type": "poly", "length": 1})
            elif isinstance(geometry, MultiPolygon):
                subpolys = []
                for polygon in geometry:
                    exterior = polygon.exterior.coords.xy
                    size = len(exterior[0]) * 2
                    i = read_count["polygons"]
                    buffers["polygons"]["coords"][
                        slice(i, i + size, 2)
                    ] = exterior[0]
                    buffers["polygons"]["coords"][
                        slice(i + 1, i + size, 2)
                    ] = exterior[1]
                    read_count["polygons"] = read_count["polygons"] + size
                    interiors = polygon.interiors
                    for interior in interiors:
                        interior_coords = interior.coords.xy
                        size = len(interior_coords[0]) * 2
                        i = read_count["polygons"]
                        buffers["polygons"]["coords"][
                            slice(i, i + size, 2)
                        ] = interior_coords[0]
                        buffers["polygons"]["coords"][
                            slice(i + 1, i + size, 2)
                        ] = interior_coords[1]
                        read_count["polygons"] = read_count["polygons"] + size
                    subpolys.append({"type": "poly", "length": 1})
                input_types.append("mpoly")
                input_lengths.append(len(geometry))
                inputs.append(
                    {
                        "type": "ml",
                        "length": len(geometry),
                        "children": subpolys,
                    }
                )
            else:
                raise NotImplementedError
        return {
            "buffers": buffers,
            "input_types": input_types,
            "input_lengths": input_lengths,
            "inputs": inputs,
        }

    def get_geoarrow_host_buffers(self) -> dict:
        """
        Returns a set of host buffers containing the geopandas object converted
        to GeoArrow format.
        """
        points_xy = []
        mpoints_xy = []
        mpoints_offsets = []
        lines_xy = []
        lines_offsets = []
        mlines = []
        polygons_xy = []
        polygons_polygons = []
        polygons_rings = []
        mpolygons = []
        buffers = self.buffers["buffers"]
        points_xy = buffers["points"]
        mpoints_xy = buffers["multipoints"]
        mpoints_offsets = self.offsets["multipoints"]
        lines_xy = buffers["lines"]
        lines_offsets = self.offsets["lines"]
        mlines = self.offsets["mlines"]
        polygons_xy = buffers["polygons"]["coords"]
        polygons_polygons = self.offsets["polygons"]["polygons"]
        polygons_rings = self.offsets["polygons"]["rings"]
        mpolygons = self.offsets["polygons"]["mpolys"]
        return {
            "points_xy": points_xy,
            "mpoints_xy": mpoints_xy,
            "mpoints_offsets": mpoints_offsets,
            "lines_xy": lines_xy,
            "lines_offsets": lines_offsets,
            "mlines": mlines,
            "polygons_xy": polygons_xy,
            "polygons_polygons": polygons_polygons,
            "polygons_rings": polygons_rings,
            "mpolygons": mpolygons,
        }

    def get_geopandas_meta(self) -> dict:
        """
        Returns the metadata that was created converting the GeoSeries into
        GeoArrow format. The metadata essentially contains the object order
        in the GeoSeries format. GeoArrow doesn't support custom orderings,
        every GeoArrow data store contains points, multipoints, lines, and
        polygons in an arbitrary order.
        """
        buffers = self.buffers
        return {
            "input_types": buffers["input_types"],
            "input_lengths": buffers["input_lengths"],
            "inputs": buffers["inputs"],
        }

1