How to convert a table to a geotable of linestrings?

[araichev]

Hi folks, Ibis newb here.
I have the attached CSV file of shape data from a GTFS feed and am trying to convert it from an Ibis table into an Ibis geotable (for lack of a better word) of linestring geometries as follows and failing.

import ibis as ib
import shapely.geometry as sg

def geometrize_shapes(shapes: ib.Table) -> ib.Table:
    """
    Given a GTFS shapes Ibis table, convert it to a geotable of LineStrings
    and return the result, which will no longer have the columns
    ``'shape_pt_sequence'``, ``'shape_pt_lon'``,
    ``'shape_pt_lat'``, and ``'shape_dist_traveled'``.
    """
    return (
        shapes
        .order_by("shape_id", "shape_pt_sequence")
        .group_by("shape_id")
        .aggregate(
            geometry=ib.literal(sg.LineString(
                zip(_.shape_pt_lon.collect().execute(), _.shape_pt_lat.collect().execute())
            ))
        )
    )

shapes = ib.read_csv(<path_to_shapes_file>)
geometrize_shapes(shapes)

This raise the following error in a Jupyter notebook.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[150], line 1
----> 1 geometrize_shapes(feed.shapes)

Cell In[149], line 16, in geometrize_shapes(shapes)
      3 def geometrize_shapes(shapes: ib.Table) -> ib.Table:
      4     """
      5     Given a GTFS shapes Ibis table, convert it to a geotable of LineStrings
      6     and return the result, which will no longer have the columns
      7     ``'shape_pt_sequence'``, ``'shape_pt_lon'``,
      8     ``'shape_pt_lat'``, and ``'shape_dist_traveled'``.
      9     """
     10     g = (
     11         shapes
     12         .order_by("shape_id", "shape_pt_sequence")
     13         .group_by("shape_id")
     14         .aggregate(
     15             geometry=ib.literal(sg.LineString(
---> 16                 zip(_.shape_pt_lon.collect().execute(), _.shape_pt_lat.collect().execute())
     17             ))
     18         )
     19     )
     21     # Optional UTM transformation
     22     if use_utm:

File ~/mrcagney/projects/gtfs_kit_next/.venv/lib/python3.12/site-packages/ibis/common/deferred.py:100, in Deferred.__iter__(self)
     99 def __iter__(self):
--> 100     raise TypeError(f"{self.__class__.__name__!r} object is not iterable")

TypeError: 'Deferred' object is not iterable

shapes.csv

How to get this right?

[gforsyth]

Hey @araichev !

So, huge caveat here that I, uh, don't really know anything about geospatial stuff, so I'd definitely check the results I've got here to make sure that they make sense.

but here goes!

First, we import a few things and read in the CSV:

[ins] In [1]: import ibis

[ins] In [2]: from ibis import _

[ins] In [3]: con = ibis.duckdb.connect(extensions=["spatial"])  # explicit DuckDB connection with geospatial extension loaded

[ins] In [4]: t = con.read_csv("/home/gil/Downloads/shapes.csv")

[ins] In [5]: ibis.options.interactive = True

[ins] In [6]: t
Out[6]: 
┏━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┓
┃ shape_id ┃ shape_pt_lat ┃ shape_pt_lon ┃ shape_pt_sequence ┃
┡━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━┩
│ string   │ float64      │ float64      │ int64             │
├──────────┼──────────────┼──────────────┼───────────────────┤
│ 1100015  │   -16.743632 │   145.668255 │             10001 │
│ 1100015  │   -16.743522 │   145.668394 │             10002 │
│ 1100015  │   -16.741258 │   145.671277 │             10003 │
│ 1100015  │   -16.744007 │   145.671074 │             10004 │
│ 1100015  │   -16.744015 │   145.671074 │             10005 │
│ 1100015  │   -16.744015 │   145.671074 │             20001 │
│ 1100015  │   -16.744732 │   145.671060 │             20002 │
│ 1100015  │   -16.747083 │   145.670999 │             20003 │
│ 1100015  │   -16.749574 │   145.671118 │             20004 │
│ 1100015  │   -16.749645 │   145.669650 │             20005 │
│ …        │            … │            … │                 … │
└──────────┴──────────────┴──────────────┴───────────────────┘

Next, we want to convert those lat/long to points instead (now here's where I might be doing this in the wrong order, but you can reverse as needed):

[ins] In [9]: t = t.select(t.shape_id, t.shape_pt_sequence, point=t.shape_pt_lat.point(t.shape_pt_lon)
         ...: )

[ins] In [10]: t
Out[10]: 
┏━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ shape_id ┃ shape_pt_sequence ┃ point                     ┃
┡━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ string   │ int64             │ point:geometry            │
├──────────┼───────────────────┼───────────────────────────┤
│ 1100015  │             10001 │ <POINT (-16.744 145.668)> │
│ 1100015  │             10002 │ <POINT (-16.744 145.668)> │
│ 1100015  │             10003 │ <POINT (-16.741 145.671)> │
│ 1100015  │             10004 │ <POINT (-16.744 145.671)> │
│ 1100015  │             10005 │ <POINT (-16.744 145.671)> │
│ 1100015  │             20001 │ <POINT (-16.744 145.671)> │
│ 1100015  │             20002 │ <POINT (-16.745 145.671)> │
│ 1100015  │             20003 │ <POINT (-16.747 145.671)> │
│ 1100015  │             20004 │ <POINT (-16.75 145.671)>  │
│ 1100015  │             20005 │ <POINT (-16.75 145.67)>   │
│ …        │                 … │ …                         │
└──────────┴───────────────────┴───────────────────────────┘

Now we want to order the columns and group by shape_id as you did above, but then pass it to the DuckDB spatial function ST_MakeLine. I couldn't find anywhere where we expose MakeLine (which seems like an oversight, but I could just be missing it), so instead we can expose it ourselves using a UDF builtin.

We need to specify the output type of the function, and make sure that the name of our "blank" function matches the DuckDB function we want to use:

[ins] In [11]: import ibis.expr.datatypes as dt

[nav] In [12]: @ibis.udf.scalar.builtin
          ...: def ST_MakeLine(*col0) -> dt.LineString:
          ...:     ...

Now we can do that order_by, group_by, and agg, and pass the collected points to MakeLine:

[nav] In [13]: t.order_by(t.shape_id, t.shape_pt_sequence).group_by(t.shape_id).aggregate(geometry=ST_MakeLi
          ...: ne(_.point.collect()))
Out[13]: 
┏━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ shape_id ┃ geometry                                                                         ┃
┡━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ string   │ linestring:geometry                                                              │
├──────────┼──────────────────────────────────────────────────────────────────────────────────┤
│ 1400001  │ <LINESTRING (-17.034 145.74, -17.033 145.74, -17.033 145.74, -17.033 145.74,...> │
│ 1420013  │ <LINESTRING (-17.026 145.728, -17.026 145.73, -17.025 145.731, -17.025 145.7...> │
│ 143W0012 │ <LINESTRING (-16.921 145.779, -16.921 145.779, -16.92 145.778, -16.92 145.77...> │
│ 150E0007 │ <LINESTRING (-17.087 145.767, -17.088 145.767, -17.088 145.766, -17.089 145....> │
│ 150E0009 │ <LINESTRING (-16.921 145.779, -16.921 145.779, -16.92 145.778, -16.92 145.77...> │
│ 1100015  │ <LINESTRING (-16.744 145.668, -16.744 145.668, -16.741 145.671, -16.744 145....> │
│ 1100023  │ <LINESTRING (-16.746 145.665, -16.744 145.667, -16.744 145.668, -16.744 145....> │
│ 110N0010 │ <LINESTRING (-16.746 145.665, -16.744 145.667, -16.744 145.668, -16.744 145....> │
│ 1110015  │ <LINESTRING (-16.791 145.681, -16.79 145.68, -16.789 145.68, -16.789 145.679...> │
│ 1110016  │ <LINESTRING (-16.92 145.779, -16.92 145.778, -16.92 145.778, -16.92 145.778,...> │
│ …        │ …                                                                                │
└──────────┴──────────────────────────────────────────────────────────────────────────────────┘

I hope that helps!

[araichev]

Thanks heaps, @gforsyth !
That works, and i've put the complete function together below for future readers of this thread.

Follow up question for you: If MakeLine were exposed in Ibis, how should i change my function?

def geometrize_shapes(shapes: ib.Table) -> ib.Table:
    """
    Given a GTFS shapes Ibis table, convert it to a geotable of LineStrings
    and return the result, which will no longer have the columns
    ``'shape_pt_sequence'``, ``'shape_pt_lon'``,
    ``'shape_pt_lat'``, and ``'shape_dist_traveled'``.
    """
    import ibis.expr.datatypes as dt

    @ibis.udf.scalar.builtin
    def ST_MakeLine(*col0) -> dt.LineString:
        ...
        
    return (
        shapes
        # Make Points as the building blocks of LineStrings
        .mutate(point=_.shape_pt_lon.point(_.shape_pt_lat))
        .order_by("shape_id", "shape_pt_sequence")
        .group_by("shape_id")
        .aggregate(geometry=ST_MakeLine(_.point.collect()))
    )

[araichev]

P.S. The the coordinate reference system, "EPSG:4326" in this case, doesn't seem to be stored in a Ibis geotable. Is that correct and therefore different from GeoPandas?

[gforsyth]

Yes, the SRID (?) isn't stored on the table by default, although there is a set_srid method. This, however, doesn't work on DuckDB as their geospatial extension doesn't support it.

[araichev]

Ah yes, the SRID. OK, good to know.

[araichev]

P.P.S. Bonus question: how to write the inverse function?
Here's an incorrect start i developed with the help of ChatGPT.
Would be quite helpful if DuckDB had an explode() method like GeoPandas does.

def ungeometrize_shapes(shapes_g: ibis.Table) -> ibis.Table:
    """
    The inverse of :func:`geometrize_shapes`.
    """
    import ibis.expr.datatypes as idt

    @ibis.udf.scalar.builtin
    def ST_AsGeoJSON(geometry: idt.Geometry) -> idt.String:
        """Converts a geometry into a GeoJSON string."""
        ...

    @ibis.udf.scalar.builtin
    def json_extract(json_str: idt.String, path: idt.String) -> idt.Array(idt.Float64):
        """Extracts an array of values from a JSON string at the given path."""
        ...

    return (
        shapes_g
        .mutate(geojson=ST_AsGeoJSON(shapes_g.geometry))
        .mutate( coordinates=json_extract(_.geojson, "$.coordinates"))
        .unnest("coordinates")   # this method does not exist in Ibis
        .select(
            shape_id=_.shape_id,
            shape_pt_sequence=_.coordinates.index.cast("int"), 
            shape_pt_lon=_.coordinates[0].cast("float64"), 
            shape_pt_lat=_.coordinates[1].cast("float64"), 
        )
    )

[gforsyth]

Wow, this was hard. I agree, an explode method would be very nice. And there might be a simpler way to do this, but for now:

import ibis
import ibis.expr.datatypes as dt
from ibis import _

@ibis.udf.scalar.builtin
def ST_AsGeoJSON(geometry) -> dt.JSON: ...  #  note return type is JSON


# Extract geo column as geo json
t = t.mutate(geojson=ST_AsGeoJSON(t.geometry))

t = (
    # Do a get of the coordinates and unwrap it as an array of float pairs
    t.mutate(
        coordinates=t.geojson.map.get("coordinates").unwrap_as("array<array<float>>")
    )
    .select(_.shape_id, _.coordinates)
    # Unnest the coordinate array with a counting offset
    .unnest("coordinates", offset="shape_pt_sequence")
    # Finally split out the coordinate lat/lon pairs into separate columns
    .select("shape_id", "shape_pt_sequence", lat=_.coordinates[0], lon=_.coordinates[1])
)

and you end up with something like

>>> t
┏━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━┓
┃ shape_id ┃ shape_pt_sequence ┃ lat        ┃ lon        ┃
┡━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━┩
│ string   │ int64             │ float64    │ float64    │
├──────────┼───────────────────┼────────────┼────────────┤
│ 1120012  │                 0 │ -16.817923 │ 145.720184 │
│ 1120012  │                 1 │ -16.817422 │ 145.720242 │
│ 1120012  │                 2 │ -16.816579 │ 145.720320 │
│ 1120012  │                 3 │ -16.815629 │ 145.720437 │
│ 1120012  │                 4 │ -16.815354 │ 145.720463 │
│ 1120012  │                 5 │ -16.813839 │ 145.720647 │
│ 1120012  │                 6 │ -16.813605 │ 145.720650 │
│ 1120012  │                 7 │ -16.812262 │ 145.720813 │
│ 1120012  │                 8 │ -16.812262 │ 145.720813 │
│ 1120012  │                 9 │ -16.811855 │ 145.720861 │
│ …        │                 … │          … │          … │
└──────────┴───────────────────┴────────────┴────────────┘

[araichev]

Nice work, @gforsyth !
Thanks for digging into Ibis's geographic dirt for me.
The current API seems tricky to use in some respects, so i might wait till it gets simpler before moving my workflows away from GeoPandas.

[araichev]

Here's the function (with lon-lat order corrected) for future readers.

def ungeometrize_shapes(shapes_g: ibis.Table) -> ibis.Table:
    """
    The inverse of :func:`geometrize_shapes`, converting back to WGS84 coordinates.
    """
    import ibis.expr.datatypes as idt

    @ibis.udf.scalar.builtin
    def ST_AsGeoJSON(geometry) -> idt.JSON: ...  #  note return type is JSON
    
    return (
        shapes_g
        # Extract geometry column as GeoJSON, which is a bit much when all we need
        # are the LineString coordinates
        .mutate(geojson=ST_AsGeoJSON(_["geometry"]))
        # Get coordinates and unwrap as an array of float pairs
        .mutate(
            coordinates=_["geojson"].map.get("coordinates").unwrap_as("array<array<float>>")
        )
        .select("shape_id", "coordinates")
        # Unnest the coordinate array with a counting offset
        .unnest("coordinates", offset="shape_pt_sequence")
        # Finally split out the coordinate lat/lon pairs into separate columns
        .select(
            "shape_id", 
            "shape_pt_sequence", 
            shape_pt_lon=_["coordinates"][0], 
            shape_pt_lat=_["coordinates"][1],
        )
    )