The earth is no longer flat - introducing support for spherical geometries in Spherely and GeoPandas

Speaker 1 [00:08]

Yeah, my talk will be about

Speaker 2 [00:10]

um

Speaker 1 [00:11]

How, at least in the Geopandas ecosystem, the Earth is no longer flat. Very quickly about myself.

Speaker 2 [00:19]

Thank you.

Speaker 1 [00:20]

I'm involved in several open source packages, so I'm one of the maintainers of GeoPandas, also of Pandas. The last years, I've worked also for my main job on Apache Arrow project, mostly the PyArrow Python bindings. But since a few months, I...

Speaker 2 [00:38]

and

Speaker 1 [00:39]

and the software engineer at Fused. So very, very briefly about what Fused is. Essentially, Fused tries to make it easier to run your Python code in the clouds. You can write your Python code, your Python function with the Python stack, the packages that you're familiar with, running it in the cloud without having to manage a cloud infrastructure yourself, and then call that function from anywhere, get the results from anywhere, integrate it very easily.

Speaker 2 [01:07]

Thank you.

Speaker 1 [01:08]

within some parallel processing caching tools. The company from where it started has some geospatial focus.

Speaker 2 [01:16]

Thank you very much.

Speaker 1 [01:18]

And that's how I ended up there. So there are some nice integrations with maps and dynamic tile servers. But in a sense, the actual server itself is not any specific to geospatial data. In the end, it's just scaling your Python code in the cloud. But back to the topic.

Speaker 2 [01:41]

of

Speaker 1 [01:42]

the talk

Speaker 2 [01:44]

Thank you.

Speaker 1 [01:45]

It's about geospatial data.

Speaker 2 [01:47]

I'm

Speaker 1 [01:49]

Maybe also just to know with the audience, like who works with geospatial data? I would like to get it. OK, so everyone, I think who is familiar with or uses regularly GeoPandas? Also large part of the room.

Speaker 2 [02:08]

Thank you.

Speaker 1 [02:09]

So you are all familiar with this kind of data. You have long-distance latitudes, maybe the reverse order. Let's not get into that discussion. And one of the tools you can use in Python to work with that kind of data when you have coordinates is Geopandas. Geopandas is based on the Shapely library. But if you use that data in GeoPandas, and Shapely specifically, that actually considers your coroner's to be

Speaker 2 [02:46]

Thank you.

Speaker 1 [02:47]

point on a plane so it's a planar coordinate so it would consider the world something like this yes you see the x and y axis which represents longitude and latitudes and you have some points of course the world as we all know is not

Speaker 2 [03:04]

Uh,

Speaker 1 [03:05]

flat like this, but we have a globe which can be approximated by a sphere.

Speaker 2 [03:11]

Thank you.

Speaker 1 [03:13]

And because of this, so you can of course work with geopandas and many people do that. You can work with geographical coordinates along with latitude. coordinates in GeoPandas, and if your GeoDataFrame is properly referenced, GeoPandas knows that your data is actually using a geographical coordinate reference system. And at that point, for example, if you would do a buffer operation, but the same will happen if you do like calculated distance, it gives a warning that your data is using a geographic coordinate reference system and the results will likely be incorrect the likely is with an understatement they're just completely incorrect

Speaker 2 [04:04]

Thank you.

Speaker 1 [04:05]

The, for example, if you, so if you calculate the distance, it will just, yeah, the distance metric will be degrees something, but that, yeah, degrees, one thing is that one degree is not the same distance everywhere on Earth. It also is just not a useful metric with which you can do anything. And so GeoPandas, for many people using GeoPandas...

Speaker 2 [04:37]

Thank you.

Speaker 1 [04:39]

will then, or the way that you can currently work with Geopandas is, or at least with those distance-based functions, is with projected coordinates. And if you have a relatively limited study area, you can use a suitable projected coordinate reference system for your study area. I'm from Belgium.

Speaker 2 [05:02]

and so

Speaker 1 [05:03]

So, for example, in Belgium, we have the Belgian Lambert 72 with its specific EPSG code. So all the data in Belgium, also a lot of open data from the government, they're all in that current reference system. Everybody's happy because that's quite accurate for the context of Belgium.

Speaker 2 [05:24]

But

Speaker 1 [05:25]

But, yeah, if you can use projected coordinates, perfect. And so at that moment, your coordinates actually are projected to a plane. You can use Geopandas and Shapely as a planner and geometry engine. Perfectly fine with only very small...

Speaker 2 [05:43]

error on that.

Speaker 1 [05:47]

But there are also situations where it's harder to use particular coordinates

Speaker 2 [05:50]

and

Speaker 1 [05:53]

first thing is that your data might already be and your data might be

Speaker 2 [05:57]

in

Speaker 1 [05:58]

in geographical coordinates, latitude longitude coordinates, and you just want to avoid any conversion of your currents because the conversion might be expensive It might be difficult to actually choose Yeah, a suitable current reference system for your data as especially if you have larger data like a larger study area or more like a global

Speaker 2 [06:24]

the

Speaker 1 [06:25]

data it's very hard to have one use one coordinate reference like projected coordinate reference system that is suitable for your entire data set so here's a typical example so whenever you project like we have our globe and you project it onto plane you always have some distortion like either it's like the shape gets distorted or the area or some other things you can never have a perfect representation. This is the famous example of the Mercator projection and it shows like when it's in full that's the Mercator projection and then it sizes down to the actual size of the country. So the Mercator projection has a very large distortion in the area of polygons when you go more towards the poles. But that's a problem if you have data you want to do global analysis on your data sets it's not that easy to yeah to just use projected coordinates

Speaker 2 [07:28]

Um

Speaker 1 [07:30]

Some other problems with projecting global data set is, for example, if you have data that crossed the, or in this case covers the South Pole, you get some kind of, yeah, edge effects. A different kind of edge problem is that The example here, the Fiji Islands, so it crosses the anti-meridian, so if you consider the Earth as this specific plane, there is a part that is entirely on the left, a part that is entirely on the right. And so if you work with planar geometries, even if you might have a single island, a single polygon, you actually have to split it to be able to represent that single island. That's already not...

Speaker 2 [08:17]

So, yeah.

Speaker 1 [08:18]

Yeah, very useful that you have to split it, but if you then do, for example, a bounding box in planner space, you get a very, yeah, not very useful bounding box.

Speaker 2 [08:34]

Thank you.

Speaker 1 [08:35]

Another example of global analysis, if you want to do something distance-based, and using a buffer here, I created the area of 500 kilometers within coast, within a distance of 500 kilometers from land. Also here, the offsets, for example, here around the equator, or more toward the North Pole, that's not actually the same distance. If you do that with... planar geometries like this Shortest path of the shortest path from the sphere. So also example of distance based calculation. So here I There is a line from I think New York to Paris whatever city in the US in Europe on The yeah projected coordinates. You might think that's the shortest path between those two points. In practice on the globe the great circle and the shortest path on the globe on the sphere is something like that. So if you plot it on this kind of projection you see and the blue line here was the the curved one but the shortest one in this projection. Final example you also have, So if you have polygons with somewhat larger edges, the question is like how, so you essentially have, this is a small polygon, so you have just a triangle, you have three points. And the question is like what do you do, like how do you connect the two points for the edge of the triangle?

Speaker 2 [10:26]

and he's

Speaker 1 [10:27]

Is this the shortest path on like the great circle as I just showed or is it? Yeah, the straight line on this representation And so if you have that triangle with those three coordinates and you consider those edges to be Planner on this representation you get a triangle like that. But if you if you would reproject that triangle using a orthographic projection looking at the poles the triangle becomes something like this and so you can already see that if you actually have those three points the real triangle that you or depending on on yeah how you created the data how you want to interpret the data the real triangle might be something like this and also here so the how you interpret edges of polygons do you interpret this as straight lines on the plane or straight lines on the sphere are two very different things because if i want to check here if a certain point is contained in that triangle and you will get a very different answer based on how you interpret your edges and so that's all um yeah related to the fact that the earth is not flat although geopana shapely will consider it as something flat and you're responsible to to ensure that your approximation, your projection is fine. A better approximation is to consider the Earth as a sphere, and that's what the libraries that I will show afterwards are doing. So using your coordinates, so your longitude-latitude coordinate, consider that to be a point on a perfect sphere. of course also this perfect sphere is still an approximation and even better approximation would be an ellipsoid where yeah because the the earth is a little bit

Speaker 2 [12:32]

And

Speaker 1 [12:33]

I'm now forgetting which direction, it will be a little bit wider than it's high, so it's not a perfect sphere, ellipsoid would be better, you can use the ellipsoid behind your coordinate reference system, but of course even an ellipsoid is always an approximation because you have local variations in... in height above the sea level.

Speaker 2 [12:58]

Um...

Speaker 1 [12:59]

But, so, by using a spherical model that is already quite accurate and is a lot simpler to implement all the algorithms of, like, does a certain polygon intersect with another polygon, all those kind of algorithms are a lot simpler to implement considering a sphere compared to an ellipsoid. And so there is a...

Speaker 2 [13:27]

But, um...

Speaker 1 [13:28]

library s2 geometry it's developed by Google it's a open source library that implements spherical geometry Yeah, algorithms. And it's, yeah, it's used by several products of Google, but it's open source and it's also... having some bindings, and for example in the R world, so I don't know if people here are familiar with the R spatial ecosystem, so you have the SF package, which is a bit the R counterpart of Teopandas. They are already using

Speaker 2 [14:07]

Um, this.

Speaker 1 [14:08]

spherical geometry engine as the default for a few years if your data are using longitude and latitude.

Speaker 2 [14:18]

Thank you.

Speaker 1 [14:19]

And so based on top of this S2 geometry library, there is another S2 geography layer.

Speaker 2 [14:28]

That's it.

Speaker 1 [14:29]

that, yeah, makes it a bit nicer to use, maps it to simple features, the point, line, string, polygons, so the simple features API that we are used to working with GeoPandas or with special databases.

Speaker 2 [14:45]

Um, and then.

Speaker 1 [14:46]

And then the Spherely library, and that's the main new library that I want to show in this talk. is a Python interface to that library. So Spherely is like Shapely, but on the sphere. And because of using being on the sphere, yeah, there are lots of those edge cases.

Speaker 2 [15:13]

like

Speaker 1 [15:14]

like being on the edges of the of the planner world and it will be handled automatically for you

Speaker 2 [15:21]

Thank you.

Speaker 1 [15:22]

The main goal of Swirly was to provide a familiar interface for people using Shapely.

Speaker 2 [15:30]

So, um...

Speaker 1 [15:32]

technical details is that it's yeah, you're using pybind11 to wrap the C++ codes from s2geometry in a Python interface and we provide some Python API Several of the functions and so there is a small code example here Using Shapely. I'm creating an array of three Shapely objects

Speaker 2 [15:58]

using

Speaker 1 [15:59]

using NumPy that is an object type array, we have three line string objects, and I can for example ask With the intersects function for each element of my array does it intersect with a certain point?

Speaker 2 [16:15]

And

Speaker 1 [16:16]

And so if you want to do the same with Spherely, so I changed the import to Spherely.

Speaker 2 [16:23]

Um, do we,

Speaker 1 [16:24]

The way to create the objects is a little bit different here. We have functions instead of the class subclass constructors because in in shape it's really we just have a single class but then we have a very similar surely dot intersects function so we try to provide a similar API as shape most of the things or yeah I have to say it's a very early like we made a 0.1 release earlier this year I think in in February.

Speaker 2 [17:01]

Thank you.

Speaker 1 [17:03]

early so it's in early stage of development but at the same time many of the basic things are actually there and are working fine so you can create all the different data types there is like you can import exports to some formats like well-known text well-known binary and you are also you can easily get data in and out Basic measurements like area distance, the basic predicates, basic overlays, they are all already available. A few notable things that are missing is actually a buffer and an area union that are, I think, two of the most used ones, at least in people using GeoPandas, that are not yet available. We certainly also have other ideas to improve it, like make it faster, make better use of the special indexing system of

Speaker 2 [17:55]

is

Speaker 1 [17:56]

s2 so s2 is based on it's a bit similar like people might be familiar with like the h3 spatial indexing system with the hexagons of different sizes at different levels so s2 has something similar with with rectangles like yeah dividing the earth in each time smaller squares and you can do very efficient like querying based on that those s2 cells that's something we haven't yet exposed in the in the api um and we also want to integrate it in uh yeah the the rest of the uh ecosystem um some examples going back to the examples that i showed earlier um so remember the triangle that is covering or not covering the uh the south pole and so if i create this triangle so with the three points

Speaker 2 [18:56]

using

Speaker 1 [18:56]

using Shapely and ask, does this triangle contain the point of the the Saipol? So Shapely is the green triangle and that's how if you consider the edges as planar straight lines, that's the green triangle and it gives false, while if you use Spherely then that creates the blue triangle and And it's a fall of my slides, I see. I don't know if I can zoom out. But yeah, you can't see it, but it shows true.

Speaker 2 [19:39]

Yeah.

Speaker 1 [19:40]

On the resolution of my screen it was visible, sorry.

Speaker 2 [19:48]

Thank you.

Speaker 1 [19:50]

Another example is area calculation and I hope that my the zoom level is a little bit so if I go for example to here to Darmstadt what I'm showing now is a web UI from Fuze that essentially and I'm whoops I'm copying lots of code here so we're here in Darmstadt

Speaker 2 [20:40]

Thank you very much.

Speaker 1 [20:42]

But I wanted to show here and I will show it here. So I'm loading some data and I did something wrong here I still paste it What I wanted to show you I'm loading some buildings data but what I want to show is that I'm essentially constructing Yeah, getting some data So the building data, I'm...

Speaker 2 [21:11]

I'm

Speaker 1 [21:12]

now converting it to Spherely, going through well-known binary, and then I want to calculate the area of those buildings in a few different ways.

Speaker 2 [21:22]

and

Speaker 1 [21:23]

And one way is, as I mentioned, if you use GeoPandas, you have to protect your data. If you don't know what is the best protection system, so one way that you can do that is there is a method in GeoPandas that is called estimate the UTM zone. So the appropriate UTM zone for your data, so it will check the total bounds of your data set and take the best approximation of the UTM zone. So we can convert it to that coordinate reference system and then calculate the area using GeoPanda, so using Shapely. I'm also calculating the area using Spherely, so using spherical coordinates, and then calculating the area, so actually the geodesic area using the ellipsoid, that's a more accurate method, using the ellipsoid behind the coordinate reference system. And so the PIPORCH package actually gives you a way to calculate the area. Yeah, based on the ellipsoid. Piper only provides that for a single geometry, so I'm doing a for loop here, as you can see, so that will be a little bit smaller. And if I zoom in more. What I wanted to show here is that, so now for this area, I'm loading the data. but it's not printing yeah it's printing a bit too much I was I should move some of those prints I wanted to let you see this part so now it's so it displayed all the buildings that I'd loaded for this area somewhere in Darmstadt and for those buildings that are here in the viewport of the map it ran this code so it ran this code for that view of the map loaded the buildings from the overture data sets for this earning box and then calculated the area using those three different ways and you can see that

Speaker 2 [23:47]

And this

Speaker 1 [23:47]

based on the spherical coordinates, there is a

Speaker 2 [23:50]

Relatively recently,

Speaker 1 [23:51]

a relative difference of like 0.3 percent so that's a typical error that you will get using spherical coordinates compared to the uh like the more accurate area based on the ellipse width while in this case using the utm zone you can see that the relative error is a little bit smaller so if you're a very small area the utm zone will be more correct it can be more accurate But, of course,

Speaker 2 [24:21]

um

Speaker 1 [24:22]

If you have global data, there is not a single UTM zone anymore or even if you have data on like full of Europe So if you want to use that method of projecting, yeah, you would have start to split your data and project separate parts

Speaker 2 [24:37]

Uh, um.

Speaker 1 [24:39]

one by one, so that becomes a lot more complex to use something like this estimate UTM zone. And so that's specifically the case where using the spherical coordinates, although there is a little bit higher error for local data, for more global data that would be a lot easier to work with. In addition

Speaker 2 [25:05]

it's

Speaker 1 [25:06]

It's also, yeah, this is essentially the code as I was showing, so using the first way, estimating the UTM zone, then using Spherely, using the ellipsoids. And in addition to the spherical methods being, if you go to more global data, the easier way, it's also actually quite fast compared to the methods based on UTM. So calculating the area itself is the blue part. That is using planar geometries is very fast. But, in this case, I also added the time to actually do the conversion of your coordinates from geographical coordinates to projected coordinates, and if you take that into account, it actually becomes slower. Of course, you might only need to do that once, and then you can continue working with your projected coordinates, but it's something to take into account if you would do that Projection just to calculate the area or calculate the distance So that's a small example of the area calculation Just to finish So what I showed up to now is all released. It's purely there is a first release you can try it out I also worked on some integration with your pandas, but that is not yet released I will just make very quickly show how it's might look like. So I'm importing a few libraries. I'm getting some land data, so a bunch of polygons.

Speaker 2 [26:51]

I'm

Speaker 1 [26:51]

I'm converting it here to, we still have to decide a bit on the exact API, whether to call it spherical geometry or call it geography. We're not yet entirely sure which way to go. But here, so I'm casting to spherical geometries. And it still looks the same, because the well-known text representation here is the same. but now you can see that the d-type is geometry spherical and if we calculate for example the area with the original data you get the warning about likely incorrect if I do that with now the spherical data we don't get a warning and we get actual value so you see that there are a lot bigger than those values in square degrees. And so now you actually get the correct area, same for if you would calculate the distance or center it, those things all work. So the idea on long term is that

Speaker 2 [28:04]

Yeah.

Speaker 1 [28:05]

Yeah, GeoPandas provides both ways, nicely integrated, all the methods stay the same, but you can say, at least initially it will be an opt-in, you can say I want to cast my geometry column to use a spherical engine under the hood, and then everything, at least everything that is implemented by Spherely will work. If you have things that don't work, it will give an error, not implement error.

Speaker 2 [28:32]

and

Speaker 1 [28:33]

Shapely has lots of methods, so Spherely will never implement everything. Buffer, though, is an example that we certainly want to have in Spherely as well. So that's, yeah, very briefly an example of how it integrates with GeoPandas.

Speaker 2 [28:51]

Yeah.

Speaker 1 [28:52]

Yeah, as I mentioned, that's a work in progress, only in a pro request for now.

Speaker 2 [28:58]

Um, yeah.

Speaker 1 [29:00]

Yeah, if you want to try it out, give it a go. So the integration with GeoPandas is not yet available, but you can use it directly already.

Speaker 2 [29:08]

Thank you.

Speaker 1 [29:10]

You can, if you load your data with GeoPandas, you can convert the Shapely geometries to Spherely, for example by going through well-known binary, and then you can get started trying it out. Yeah, a few people to thank. So I, together with Benoit, we worked on the Spherely package, also thanks to a NumFocus small development grant.

Speaker 2 [29:35]

Thank you.

Speaker 1 [29:36]

And thanks to the other people involved, and also certainly the jury and Edzer from our special community, who actually started doing this a few years ago. And we are now being able to benefit from their work to bring that to Python as well. Okay, so that was my presentation. You can find, if you want to see the slides, there is a link here to the GitHub repo with the slides.

Speaker 3 [30:20]

Wow, thank you. Jerry's that was a great talk. I pretty much enjoyed everything and I believe everybody did So let's go on a slider to Look at some questions so the first question is always from anonymous is Do you plan to support changes from a sphere to a squeeze sphere and allowing for a 3d dimension like altitude?

Speaker 2 [30:49]

Thank you.

Speaker 1 [30:50]

The first part I didn't entirely understand about the squeezed geometry. But so I wanted to answer the second part. So altitudes, that's something that is not supported by the S2 geometry library. They only do essentially on the unit vector on the globe. And yeah, so it's only longitude latitudes, no altitudes. So there we are limited by what the underlying library does.

Speaker 2 [31:17]

And for.

Speaker 1 [31:18]

For example, Shapely does support some form of altitude like with z and also m coordinates But it's also very limited. Many of the algorithms in Shapely will only consider actually x and y and not your Altitude so that's in general in like geopandas quite quite limited

Speaker 3 [31:43]

The second question is, does MCP support servers calling the LLM in reaction to an event rather than only the LLM calling the server?

Speaker 1 [31:53]

I suppose that's someone who put it in the wrong room.

Speaker 3 [31:53]

I suppose. Okay, we have two more. Is there a chance that the methods of Sferi will be added to Shapely and therefore to your pandas?

Speaker 2 [32:09]

Um, so...

Speaker 1 [32:10]

So we don't plan to add it, actually, to Shapely. Shapely is really based on the underlying Geos library. That's the underlying C++ library that many packages use. Also, if you're using Post.js, it's using Geos. If you're using R, if you're using QGIS, they're all using the Geos library for the actual operation. So Shapely is limited to Geos, and Geos is Planner. That's their scope. and that's the reason that we created this separate package spherely to focus on spherical geometry engine

Speaker 2 [32:45]

Thank you.

Speaker 1 [32:46]

So Shapely and Spherely will always be separate and if you want to use it directly you will have to do that consciously and Use one or the other and you can convert from one to the other by going for example through well-known binary We do want to integrate it as I showed in the end. We want to integrate it in GeoPandas that at least for GeoPandas users Where you might not use Shapely that much directly

Speaker 2 [33:11]

uh, in

Speaker 1 [33:12]

becomes a quite fluent experience like you have your geometry column and you can say like now convert it to spherical like use a spherical engine instead of the planar engine and that that becomes quite yeah seamless integration but that means that if you interact with single geometries of your two data frame like you get one polygon out of it you will get this really polygon instead of shapely polygon and that will yeah behave slightly differently like they are not exactly the same

Speaker 3 [33:46]

Okay, we have two more questions is the API of very similar to that of shapely can you can Can use a stropping replacement does it also use same dependencies like Geo's library under the hood?

Speaker 1 [34:02]

Yeah, so I partly mentioned that, I think, in answering the last question. So it's not using geos, it's using a different neighborhood library. But it's very similar in that it's also, at least, it integrates very well with NumPy. So if you work with an array of geometries. So here I had an example, and you see that it's not exactly a drop-in replacement. Because we made some decisions in Spherely to not have these subclasses. So Shapely has a point, line string, polygon subclasses, and you can use them as constructors. In SphereD, we didn't implement all the subclasses, so we have a createLineString function instead. That's a difference in this example, but then once you have your geometries created or you converted them from well-known text, and at that point, it's very similar. Then here, the function intersects works exactly the same. They both work with single geometries or with arrays of geometries.

Speaker 2 [35:07]

Thank you.

Speaker 1 [35:08]

So not exactly drop-in, but we want to make it at least very familiar for people that know Shapely to use Spherely.

Speaker 3 [35:20]

Great. The last question, no? We have, I think we can, can we take two more questions? Okay. So the surface of the Earth is still flat in spherae, like Belgium and Ennard, actually.

Speaker 1 [35:35]

Sorry, can you say it again?

Speaker 3 [35:39]

The person is just asking if the surface of the earth is still flat in theory like in Belgium or Netherlands, I don't know if this one is the series of somebody's just joking

Speaker 1 [35:50]

Yeah, it's true.

Speaker 3 [35:51]

Anybody ask this question, can you expand more?

Speaker 1 [35:56]

It's still a perfect sphere. So yeah, everything is locally everything is flat

Speaker 2 [36:00]

That's it.

Speaker 1 [36:01]

That's true.

Speaker 3 [36:02]

Okay. The last question is from Anonymous. He or she is saying, great work. Thank you.

Speaker 1 [36:11]

Thank you.