Converting lat/lon numpy array(s) to geoarrow-rust objects?

[wx4stg]

Hello! I have a notebook that reads in weather radar data, does some coordinate transformation, and produces a numpy array of shape (366143, 4, 2), representing 300k polygons that have four corners (they're not quite square...) and each corner has latitude and longitude.

To plot this data with lonboard, I'm creating POLYGON (( ... )) strings, (using numpy string operations rather than shapely polygons to avoid the overhead of a python for loop), then using geoarrow-pyarrow to convert this array of strings to a geoarrow object, then converting to a geodataframe and finally plotting.

The relevant lines of code below are the monstrosity that are my numpy string operations and the calls to geoarrow.pyarrow.as_geoarrow and geoarrow.pyarrow.to_geopandas...

Minimal example

import numpy as np
import geoarrow.pyarrow as ga
import geopandas as gpd

# poly_coords is an array of shape (366143, 4, 2)

lon_lat_pairs = np.strings.add(np.strings.add(poly_coords[:, :, 0].astype(str), ' '), poly_coords[:, :, 1].astype(str))
polygon_strings = np.strings.add(np.strings.add(
    np.strings.add(np.strings.add(
        np.strings.add(np.strings.add(
            np.strings.add(np.strings.add(
                lon_lat_pairs[:, 0], ', '), lon_lat_pairs[:, 1]), ', '), lon_lat_pairs[:, 2]), ', '), lon_lat_pairs[:, 3]), ', '), lon_lat_pairs[:, 0])

polygon_strings = np.char.add(np.char.add('POLYGON ((', polygon_strings), '))')
geom_arrow = ga.as_geoarrow(polygon_strings)
geometry = gpd.GeoDataFrame(geometry=ga.to_geopandas(geom_arrow))
# Then I plot it as a SolidPolygonLayer... you get the idea

This works "great"!

However, I'd like to use this with upcoming versions of pyodide, geoarrow.pyarrow is not supported, but geoarrow-rust is!

---

Is there a way to replace my calls to gearrow.pyarrow with geoarrow-rust to create a geopandas dataframe with the geometry of these points (quickly, without a python for loop)?

And, I'm reaching here, is there a way to do this using just the original numpy array, without the intermediate WKT string representation?

[kylebarron]

First, super cool to see weather (radar?) data visualized in this way:

If you have an open, public-facing example, we should add it to our collection of user examples in the docs website.

Otherwise, I have a few comments:

To plot this data with lonboard, I'm creating POLYGON (( ... )) strings

😬 that sounds like a mess 😅

using numpy string operations rather than shapely polygons to avoid the overhead of a python for loop

As an aside, note that even with shapely you don't need a Python for loop here. You can use the shapely.polygons constructor with your numpy array input directly.

I will admit that that constructor isn't the easiest to use. It always takes me some time to recall how to use it, but it is possible to use with direct numpy input.

Still though, that gives you vectorized creation of shapely geometries, and lonboard would still need to convert those geometries to GeoArrow, so a GeoArrow-native approach would be faster.

However, I'd like to use this with upcoming versions of pyodide, geoarrow.pyarrow is not supported, but geoarrow-rust is!

To be clear, geoarrow-pyarrow can still be used in any normal Python environment outside of the browser. And pyarrow actually now works in pyarrow, the missing piece is just the geoarrow-c dependency that would need to be built for pyodide. I imagine that's not too difficult but you'd have to ask Dewey on that repo.

Is there a way to replace my calls to gearrow.pyarrow with geoarrow-rust to create a geopandas dataframe with the geometry of these points (quickly, without a python for loop)?

Well, with geoarrow-rust you need to work with Arrow tables, not GeoPandas GeoDataFrames. But in principle yes. The smallest change from your existing workflow is to create an Arrow string column from your numpy string data (e.g. with pyarrow.array(numpy_string_array)) and then pass it to from_wkt, which will parse it to GeoArrow in a vectorized way.

Dewey created geoarrow-pandas and geoarrow-pyarrow; geoarrow-pandas currently only works with geoarrow-pyarrow, not geoarrow-rust, so you can't use a GeoPandas GeoDataFrame with a geoarrow-rust column. It would be cool to implement support for geoarrow-rust as well, but I don't have time for that right now.

And, I'm reaching here, is there a way to do this using just the original numpy array, without the intermediate WKT string representation?

Yeah you really don't want to go through WKT. I haven't implemented native constructors in geoarrow-rs like the shapely.polygons function because I'm not sure what is the cleanest API to present (while still being vectorized).

You can still do this today in an "advanced" way with a bit of knowledge of the GeoArrow format. Assuming that you have a numpy array of shape (N, 4, 2), then every polygon is "simple", which makes this easier. Note that polygons must be closed so copy the [:, 0, :] axis into [:, 4, :] to get an array of shape (N, 5, 2).

Then reshape your data to a flat array of all coordinates ordered as XYXYXY. I'm 99% sure all you need is an array.ravel("C").

Then construct an offsets buffer for each ring. Since every ring has exactly 5 coordinates, this is just np.arange(0, (N + 1) * 5, 5, dtype=np.int32). You want an array of [0, 5, 10, 15, etc]. This is the offset of each ring, so the first ring has 5 coordinates, and so on. You need N + 1 because the offsets buffer should be one longer than the number of rings.

Then construct an offsets buffer for each geometry. Since every polygon has exactly 1 ring, this is just np.arange(N + 1, dtype=np.int32).

Then you can construct a geoarrow array as such:

lonboard/lonboard/_geoarrow/extension_types.py

Lines 225 to 237 in 7562834

	arrow_coords = fixed_size_list_array(
	Array.from_numpy(coords.ravel("C")), len(dims)
	)
	arrow_rings = list_array(ring_offsets, arrow_coords)
	arrow_geoms = list_array(geom_offsets, arrow_rings)
	extension_metadata["ARROW:extension:name"] = "geoarrow.polygon"
	field = Field(
	field_name,
	arrow_geoms.type,
	nullable=True,
	metadata=extension_metadata,
	)
	return field, arrow_geoms

(You can look back in that file's history if you want to see how to do it with pyarrow)

I hope this is somewhat easy to follow!

[wx4stg]

First, super cool to see weather (radar?) data visualized in this way
If you have an open, public-facing example, we should add it to our collection of user examples in the docs website.

It is indeed weather radar data, using openradar/xradar. @mgrover1 (one of the authors of that library) and I attended your SciPy talk and both immediately realized the fantastic applications of using lonboard for meteorological data viz.

I will absolutely be contributing an example here (and maybe there too with some more science-focused examples), but I'd like to polish things a little bit (for obvious reasons, I don't particularly want my example unnecessarily doing something as silly as "putting all the coordinates of radar gates into strings", haha). If I get more comfortable/proficient with lonboard, there's also a repository of notebooks with more in-depth trainings for geoscientists, Project Pythia, and I think they'd appreciate lonboard as well. This is a really awesome project and y'alls efforts are very much useful in a lot of fields.

---

I think you'd be surprised at how fast the string operations run, including I/O of the radar data it takes 6 seconds to generate a plot on a pretty modest laptop (intel 1240P, integrated graphics only). But I do agree that going through WKT is extremely silly and unnecessary overhead.

You can still do this today in an "advanced" way with a bit of knowledge of the GeoArrow format.

I will admit this is essentially my first experience with geoarrow. My understanding is that Arrow/Parquet as being an efficient/binary replacement for the csv, and the logical extension is that adding a "geometry" column (either WKT or WKB) to a parquet... then you have a "geoarrow".

# Copy the first coordinate point to the end of the array
poly_coords_closed = np.concat((poly_coords, poly_coords[:, 0:1, :]), axis=1)
# Shape of poly_coords_closed is (366143, 5, 2)
# Shape of np.ravel(poly_coords_closed) is (3661430,)
poly_coords_arrow = fixed_size_list_array(Array.from_numpy(np.ravel(poly_coords_closed)), 2) # I'm a bit unsure if this should be 2 or 5 or 10, see below
ring_offsets = Array.from_numpy(np.arange(0, (poly_coords_closed.shape[0] + 1), 5, dtype=np.int32))
arrow_rings = list_array(ring_offsets, poly_coords_arrow)
geom_offsets = Array.from_numpy(np.arange((poly_coords_closed.shape[0]/5), dtype=np.int32)) # I had to divide this by 5, which really feels wrong...
arrow_geoms = list_array(geom_offsets, arrow_rings)

extension_metadata = {
    'ARROW:extension:name' : 'geoarrow.polygon'
}
field = Field('geometry', arrow_geoms.type, nullable=True, metadata=extension_metadata)
layer = SolidPolygonLayer(table=arrow_geoms)
m = Map(layer, basemap_style=CartoBasemap.Positron)

So, I'm not sure whether the poly_coords_arrow list_size should be 2 or 5 (or 10?). one list per value (2), one list per lat/lon pair per polygon (5), or one list for each polygon (10 values total?)

I also had to go against your suggestion of using an arange from 0->N+1 for the geom_offsets, it gave an error of Exception: Invalid argument error: Max offset of 366143 exceeds length of values 73228 when creating the arrow_geoms array.

However, neither arrow_geoms nor field have a schema attribute, so passing them to the "table" kwarg of a solidpolygonlayer doesn't work.

I've uploaded my notebook if you're curious: https://gist.github.com/wx4stg/f7f05acfe7e6ecafceb5738164e06f81

[kylebarron]

It is indeed weather radar data, using openradar/xradar. @mgrover1 (one of the authors of that library) and I attended your SciPy talk and both immediately realized the fantastic applications of using lonboard for meteorological data viz.

🙌

I will absolutely be contributing an example here

Just as a quick note, the current thinking is that for notebooks hosted in the lonboard repo we want easy-to-follow and very general examples, but for more complex notebooks to have links in the docs to other people's externally-hosted notebooks. See https://developmentseed.org/lonboard/latest/examples/#third-party-showcase. That lessens the maintenance burden for "Lonboard-official" examples, especially for stuff like radar data that might need domain-specific knowledge that I don't have.

I think you'd be surprised at how fast the string operations run

It's impressive that it's that fast (really an example of numpy v2's great string support) but apart from that in my experience it's easy to create malformed WKT. I'd definitely like to figure out the right answer that is both easy for you/users and doesn't go through WKT.

I will admit this is essentially my first experience with geoarrow.

Just to be clear, I'm mentioning this as a way to unblock you now, and also in case you're curious in learning a bit more how GeoArrow works under the hood, but I'd also like to work towards simpler ways to construct GeoArrow data from raw coordinates.

My understanding is that Arrow/Parquet as being an efficient/binary replacement for the csv, and the logical extension is that adding a "geometry" column (either WKT or WKB) to a parquet... then you have a "geoarrow".

Yes... but the added important distinction is that there are two "kinds" of GeoArrow "native" and "serialized" encodings. WKT and WKB are serialized encodings, which means that the data still has to be parsed to be used. Lonboard accepts serialized encodings, but parses it to the GeoArrow native encoding, which is sent to the browser. You can tag an Arrow string array as a geoarrow.wkt array, and that will work out of the box in lonboard, but with some overhead because lonboard will currently use shapely to parse that WKT and then create GeoArrow from it.

The instructions I tried to give you earlier were for creating a "native" PolygonArray. In the native arrays, there's O(1) access to any coordinate, so it's as efficient as possible.

Let me take a look at your notebook

[kylebarron]

ring_offsets = Array.from_numpy(np.arange(0, (poly_coords_closed.shape[0] + 1), 5, dtype=np.int32))

You missed the * 5 in the second arg, so you're not creating the right length of offsets array.

You also need to pass in a table, not just a geometry to the layer class, because most of the time you have some attribute per geometry. This is something maybe we could relax in the future.

See repro notebook here: https://github.com/kylebarron/lonboard-radar-example/blob/d372b6526d6fd766e6767c80ccbfea27914590b2/radar-viz-lonboard.ipynb. That repo also has a pixi env file if you want to reproduce it.

[wx4stg]

Ah, yes that's what I was missing. Thanks for looking so far deeply into this! Excited for what can be made with lonboard

[wx4stg]

continuing discussion for weather radar viz at ARM-DOE/pyart#1653