Introducing first-class Spatial Data Types in Apache Parquet

import geopandas


# ~17s
df = geopandas.read_file("ns-water_elevation.fgb")
df.info()
#> <class 'geopandas.geodataframe.GeoDataFrame'>
#> RangeIndex: 11438004 entries, 0 to 11438003
#> Data columns (total 4 columns):
#> ...
#> dtypes: float64(1), geometry(1), object(2)

import pyarrow as pa
from pyarrow import parquet

# ~0.5 s
table = parquet.read_table("ns-water_elevation.parquet")

geopandas.GeoDataFrame(table.to_pandas()).geometry
#> AttributeError: There are no existing columns with geometry data type.
#> You can add a geometry column as the active geometry column with df.set_geometry.

import pandas as pd

df["geometry"] = df.geometry.to_wkb()
#> UserWarning: Geometry column does not contain geometry.
df_pd = pd.DataFrame(df)
df_pd.to_parquet("foofy.parqeuet")

df2 = table.to_pandas()
df.geometry.geom_equals(df2.geometry)
#> TypeError: (<class 'geopandas.geoseries.GeoSeries'>, <class 'pandas.core.series.Series'>)

df2 = table.to_pandas()
df.geometry.geom_equals(geopandas.GeoSeries.from_wkb(df2.geometry))
#> UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
#> Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.

f = parquet.ParquetFile("ns-water_elevation_geo.parquet")
f.schema_arrow.metadata
#> {
#>   'geo': '{
#>      "columns": {
#>        "geometry": {
#>          "encoding": "WKB",
#>          "crs": {...long fancy but very precise value...},
#>          "bbox": [228907, 4807442, 765782, 5234427]
#>        }
#>      },
#>      ...
#>   }'
#> }

from geoarrow.pyarrow import io

# ~0.5s
table = io.read_geoparquet_table("ns-water_elevation_geo.parquet")
table["geometry"].type.crs
#> ProjJsonCrs({"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "Compo)

from sedona.spark import *

config = SedonaContext.builder().getOrCreate()
sedona = SedonaContext.create(config)

df = sedona.read.format("parquet").load("s3://...")
df.sparkStuff()

ax = df[df.ZVALUE > 500].plot()

url = "https://raw.githubusercontent.com/geoarrow/geoarrow-data/v0.2.0/natural-earth/files/natural-earth_countries.fgb"
countries = geopandas.read_file(url)
countries = countries[countries.name == "Canada"].to_crs(df.crs).intersection(df.union_all().convex_hull)
ax = df[df.ZVALUE > 500].plot()
countries.boundary.plot(color="purple", ax=ax)

import geoarrow.pyarrow

f = "ns-water_elevation.parquet"
tab = parquet.read_table(f)
geopandas.GeoDataFrame.from_arrow(tab).info()
#> <class 'geopandas.geodataframe.GeoDataFrame'>
#> RangeIndex: 11438004 entries, 0 to 11438003
#> Data columns (total 4 columns):
#>  #   Column     Dtype
#> ---  ------     -----
#>  0   FEAT_CODE  object
#>  1   FEAT_DESC  object
#>  2   ZVALUE     float64
#>  3   geometry   geometry

import duckdb

duckdb.load_extension("spatial")
duckdb.sql("CALL register_geoarrow_extensions()")

tab = parquet.read_table(f)
duckdb.sql("SELECT geometry FROM tab LIMIT 5")
#> ┌──────────────────────────────────────────────────────────────────┐
#> │                             geometry                             │
#> │                             geometry                             │
#> ├──────────────────────────────────────────────────────────────────┤
#> │ POINT Z (255919.02740000002 4831906.580700001 6.699999999999363) │
#> │ POINT Z (256023.9274000004 4831899.6807 1.399999999999977)       │
#> │ POINT Z (255988.62739999965 4831942.2807 1.399999999999977)      │
#> │ POINT Z (255900.7274000002 4831946.8807 1.399999999999977)       │
#> │ POINT Z (255862.2274000002 4831922.580700001 9.900000000000091)  │
#> └──────────────────────────────────────────────────────────────────┘

sd.read_parquet(f).to_view("tab")
sd.sql("SELECT * FROM tab").to_pandas()