Skip to content

2026

V0.4.1 Updates

Highlights

  • All geometry WKB will be stored in the df as little endian ISO / Properly handling input wkb with SRID added to it 🐛
  • No longer renaming your geometry column "geometry" without your permission 🚫
  • scan_spatial will no longer work with geoparquet 😢
  • scan_geoparquet added to do that! 😄 🎉

Little Endian ISO WKB

In V0.4.0 new functionality was added to user polars expressions to parse geometry WKB and performa operations without converting to shapely first. This worked, but parsing different flavors/byte orders looked like it was going to be a lot more work and less efficient than only dealing with one flavor/byte order. Now when the WKB is added to the struct, if it is not already little endian/ISO it will be automatically transformed to it with shapley. This will also handle a bug I had when working with coordinates of extended WKB points with an embedded SRID.

No more automagically changing names of gometry columns to "geometry"

Not all geospatial sources call their geometry field "geometry". Spatial polars shouldn't have renamed the column if the data had a different name for the field. It no longer does that.

Scan spatial/Scan geoparquet

When Spatial polars started, polars did not like to read parquet files that had extensions on columns. To work around that, spatial polars used to read the batches of geoparquet data with pyarrow and then convert that to a polars series in the scan_spatial function. Polars can now handle reading in geoparquet files without issue, that said, the geoparquet metadata with the CRS of the data is lost when you read a geoparqut file with polars. I've added a new scan_geoparquet function to spatial polars that will call polars' scan_parquet but also finds all(1) the geometry columns and converts them to the spatial polars geometry struct using the crs from the geoparquet file's column metadata.

  1. The old scan_spatial function only converted the column that the metadata designated as the primary column
  • in wkb
  • 6 min read

Polars, Meet WKB

Prior to V0.4.0 spatial polars treated WKB in the geometry structs as a black box. For all spatial expressions, the WKB was shipped over to shapely, converted to a shapely geometry object, then whatever operation you requested was performed on the shapely goemetry object, then the results were passed back to polars. It was simple, and shapely's numpy ufuncs certainly aren't what I consider slow, it was all good! Then DuckDB and SedonaDB came in with their blazing speeds, and made me once again reconsider my concept of slow 🐢.

An Idea Enters 💡

A few weeks ago, I stumbled on the binary expressions documentation offered by polars. I'm sure I've seen that page before, but for some reason when I was there, it clicked! A lot of spatial operations(1) are relatively simple math problems. Example: to know the distance between two points(2), we just need to use the pythagorean theorem, which is trivial to implement using polars expressions. I had a hunch that if we were to use polars expressions to parse the binary to extract the X/Y coordinates of the points and then supply those coords to polars expressions to compute the distance, we could get the same answer as before, but hopefully quicker, because we wouldn't need to create an entire shapely geometry object just to do some simple math.

  1. Especially ones involving points
  2. In cartesian space (which is what shapely does)

Let's Do it!

Armed with the idea, I set forth to see what I could do. I took a dataframe with a bunch of points, and wrote a bare bones expression to parse the X coordinate from the WKB. It worked super fast, read WAY faster than sending the WKB to shapely and asking shapely to give me the X coordinate. I quickly spun up another expression to return the Y coordinate and then another to compute the pythagorean theorem to determine the distance to an arbitrary other point, and finally if the distance was less than an input distance to emulate the dwithin function of shapely. Once again it was overwhelmingly faster than the shapely based function, I was extatic 😳! Excitement aside, I knew what I was prototyping with was only dealing with little endian ISO WKB. The expression I ultimately wanted needed to be taking into account WKB with big endian byte order, little endian extended WKB, and big endian extended WKB(1), and would still need a fallback (at least for now) to shapely for non-point geometry types. Although I knew I needed to write a lot more code, I figured with the parallel nature of polars it would be a bit slower to take into account all those possibilities, but was hopeful that it would still net a large perfomance gain.

  1. At the time I wasn't even thinking about extended WKB with an embedded SRID... which would further complicate matters... more on that later :-/

Make it Robust

I began chaining out a big pl.when().then().otherwise() expression to handle the different byte orders and WKB flavors. I set up an expression to parse the binary to determine which geometry types were in the series, if they were little endian points, I supplied the WKB to my dwithin expression for little endian points, big endian points went to another expression, otherwise use shapely to compute the distance. Long story short, that didn't pan out. There's a big ole warning in the polars.when docs that clearly states that all the expressions would be evaluated in parallel(1), so even though the data I was experimenting with was all points, I after adding all the different .when().then() chains and having the fallback to use shapely, it was still doing all the computations including the fallback to shapely even though it wasnt going to use that result since I had the little endian point operation first in my chain... It was slower than just using shapely. I was sad, but I didnt give up!

  1. I didn't honestly pay attention to the big warning note there... I'm the problem, it's me.

Where to Now?

I deicded to use a different expression first to split the big/little endiain point WKB data into different fields of struct, and then the non-points to a third field of a struct, so each row would have one of the three fields populated leaving the other two fields null. I passed the field of little endian WKB to one expression, the big endian to another, and the non-point geometry types field would go to shapely for the computation, then coalescing the results for each field in the struct. In my experimentation with only points, since the data being passed to shapely was all None, shapely returned None back as a result super fast(1), and the little/big endian expressions were being run quickly, so overall I was getting much better performance(2) without touching the public API(3).

  1. I'm guessing it's got some sort of no-op shortcut when you give it a None to work from, it was almost instant on 6M points
  2. Not as fast as just having only one expression which only worked for the little endian points, but still loads faster than the shapely method
  3. I didnt want to make some sort of "only_points_dwithin" function or something like that, that feels complicated and confusing to me

Computative Results

I pushed the updates up to github/pypi as V0.4.0, and asked the crew over at spatial bench if they'd be willing to kick off the github action to run that benchmark suite so I could see the results there. Previously Q1 of the spatial bench benchmark, which uses the dwithin function, was taking spatial polars in the 3.8-3.9s range. Now running spatial polars v0.4.0 the exact same query, producing the same results only took 0.82s 🚀 that's still a fair bit slower than DuckDB or SedonaDB, but it's not totally getting left in the dust anymore!

Follow On

Overall I'm pretty jazzed 🎺 about the results, and think I'll mess around with adding some more polars binary based expressions in the future. The more I was thought about it though, dealing with the different byte order and WKB flavors seems like a lot of overhead. Not to mention I also woke up out of nowhere one night this week realizing that I hadn't accounted for embedded SRIDs, so results coming from extended WKB with an embedded SRID are wrong(1). In version 0.4.1 I've decided to take all WKB input to the geometry series, and if it's not already little endian ISO, it will be converted to it, which resolves that SRID bug 🐛. This will also ease the maintenance burden of all the different varieties of WKB, and shorten the when.then.otherwise chains, so theoretically we should get some more performance, provided the input data is already little endian ISO.

  1. I'm being 💯 serious... who thinks of this stuff in their sleep??? ... me I guess... 🤷
  • in join
  • 7 min read

Lazy Spatial Joins

When I initially started spatial polars, one of the things I really didn't like about what I had created was that it wasn't capable of doing a spatial join from a lazyframe. I messed around with some attempts at registering a spatial lazyframe namespace/adding spatial join functionality to it. My original attempt involved collecting the data from the lazy frames within the join function. The problem with that approach, is that it was collecting the query to build the spatial join's intermediate frame which joins the two input frames immediately, before the rest of the query was executed. That defeated the entire idea of the lazy spatial join. Clearly, I wasn't going to go with that. I knew there likely was a way to do what I wanted, but how to implement it evaded my thoughts until I bumped into the pl.defer() function a few weeks ago. Defer was added to the polars API over a year ago, I'm not sure if I missed it in the release notes when it got added, or if I simply didn't realize how useful it could/would be, but it's relatively new to me, and it's GREAT 💥!

Tell me what you're planning, I'll do it later.

Defer, per the docs, Takes a function that produces a DataFrame but defers execution until the LazyFrame is collected. It's absolutely perfect for what I needed! Inside the lazyframe's join function, I have another function wrapped with defer which collects the geometries of the left and right frames, then builds out a STRtree, and transforms the results into an intermediate dataframe that is used to join the left and right frames together. Because that joining frame is created by a function called with pl.defer, it doesn't get executed until the query is actually collected! lazy spatial join problem solved! well kinda...

What can't it do?

Streaming from traditional GIS sources

So one of the really cool things about polars lazyframes and the streaming engine is it's ability to work with datasets that are bigger than the RAM on the box processing the query. If you can't fit your data into RAM, this lazy spatial join still wont work, it needs to know all the geometries from both frames at the same time. There's no spilling to disk or anything like that. For the datasets I've worked with for the past 20+ years, that's not an issue, but if you've got super huge data you're working with, and/or limited RAM you're still out of luck here. That's not an issue I run into much/ever, so I'm honestly not motivated to try to solve it(1). I like being upfront and realistic 😄. That out of the way, I do like the idea of being able to use the polars streaming engine to write batches of data, so perhaps I'll try and figure it out at some point 🤔(2). Worth noting here, if your data is in parquet format, not a more "traditional" format (eg shapefiles, file geodatabase, geojson... stuff that spatial polars uses pyogrio to access) in my experimentation thus far, collecting queries which involve a lazy spatial join with the streaming engine works as expected.

  1. Just being real 😄, if you're in need of spatial ops on billions of rows, I doubt this will ever be the appropriate package
  2. I think there's something going sideways between polars and pyogrio...

So what CAN it do?

Regular spatial joins, join nearest and centroid KNN joins are all now available for your spatial data processing delight.

The syntax for the joins are the same as for the normal dataframes, we'll just operate on lazyframes instead of dataframes:

Lazy Spatial Join
import polars as pl
from spatial_polars import scan_spatial

lake_lf = scan_spatial("https://naciscdn.org/naturalearth/110m/physical/ne_110m_lakes.zip") # (1)!
boundary_lf = scan_spatial("https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip")  # (2)!

lake_boundary_lf = (
    lake_lf.spatial.join(  # (3)!
        other=boundary_lf,  # (4)!
        how="inner",  # (5)!
        predicate="intersects",  # (6)!
        on="geometry",  # (7)!
        suffix="_boundary",  # (8)!
    )
    .select(
        pl.col("name"),  # (9)!
        pl.col("SOVEREIGNT"),
        pl.col("geometry"),
        pl.col("geometry_boundary"),
    )
    .sort("name")  # (10)!
)

lake_boundary_df = lake_boundary_lf.collect(engine="in-memory")  # (11)!
print(lake_boundary_df)
  1. Scanning the lakes
  2. Scanning the boundaries
  3. Starting with the lake_df we'll start our spatial join
  4. Specifying to join the lakes to this boundary_df

  5. We'll use an inner join so as to only return rows for lakes that actually intersect a boundary. If a lake does not intersect a boundary polygon we won't have a row for that lake in our output dataframe. Likewise, if a boundary doesn't intersect a lake, the resulting dataframe won't have a row for that boundary.

    Note

    This could have been left off, because how='inner' is the default

  6. Use the 'intersects' spatial predicate so if any part of the lake shares any space with the boundary we'll join the lake to the boundary. Since we've specified an inner join, if a lake intersects more than one boundary, we'll get more than one row for the lake since it's joined to more than one boundary.

    Note

    This could have been left off, because predicate='intersects' is the default

  7. Since the name of the geometry struct is 'geometry' in both of our dataframes, we will use the on parameter, if we wanted to use a different column name for each of the dataframes we could use the left_on or right_on parameters.

    Note

    This could have been left off, because on='geometry' is the default

  8. Since we're joining the dataframes with a common column name (geometry), a suffix must be applied to the columns of the right dataframe that have names that exist in the left fram, because we can't have two columns with the same name. we'll use "_boundary" as the suffix to clarify that the geometry of the right frame came from the boundaries dataframe.

  9. Selecting the columns to make the lake name and SOVEREIGNT columns show up before the lake and boundary geometry columns.
  10. Sort by the lake name just to make the results look nice in our output dataframe.
  11. Note that "in-memory" engine.. streaming isn't working here for some reason. not sure what's causing it, but in my testing thus far my python kernel crashes when I try to use the streaming engine. I suspect something with pyogrio and polars aren't playing well together here, but haven't gotten to the root of the problem yet.

Why lazy?

When we execute a lazy query, we gain the optimizations from the polars query engine. In the example above, the lakes data has 39 columns, and the boundaries has 170. The result of our query only returns two string columns, and the geometries from the lakes and the boundary. If we look at the graph for the query, we can see that polars sees that we only need a two columns from each of our sources!(1), so the query engine tells our io-plugin to only read the name/geometry from the lakes, and the SOVEREIGNT/geometry from the boundaries. If we had to read all the columns from each of the source and then do our join and only keep a few fields at the end that's just wasting time/resources.

  1. Note the PYTHON SCAN 2/38 and 2/169 boxes at the bottom of the graph indicating only two columns will be read from the sources.

Lazy spatial join show graph output

That's it for now

I really like building stuff, and I enjoy writing up blogs like this about the stuff I've created. I do NOT like writing the end of them though. Like, what do I write here? I'm not trying to sell anything(1), I'm not necessarialy trying to get anyone to use this package(2), I don't have some sort of closing slogan(3). I guess I'll just channel my high school English class informative writing coursework, and tell you what I've already told you. Spatial polars can now do a spatial join from a lazyframe without the need to collect the query into a dataframe first! 🥳 (4)

  1. It's free.
  2. I do think it's useful, and hopefully you do too!
  3. Ohhhh... maybe I need some sort of closing line like when Matt O-Dowd ends episodes of spacetime with a cleverly worded sentence that ends with "Spacetime", and then the episode ends, that'd be pretty rad, but I don't think I'm clever enough to pull off a lot of different ways of plugging "spatial polars" into the end of a sentence...
  4. High school english class didnt tell me anything about using emojis in my informative essays... because... well they didn't exist yet 😆 👴