geopython / geopython/pygeofilter

support for DuckDB

Open
#90 5 comments 0 reactions 0 assignees View on GitHub
Dominant language
Python
Stars
90
Forks
44
PR merge metrics
No merged PRs in 30d

Description

Hi,
I am currently trying to figure out a way to use `pygeofilter` for converting CQL2 to a DuckDB query using the [spatial extension](https://github.com/duckdb/duckdb_spatial). This extension is quite prototypical and I stumble across a couple of caveats. Here's my current approach:
```python
import duckdb

duckdb.install_extension('spatial')
duckdb.load_extension('spatial')
```
Define the CQL2 filter
```python
cql2_filter = {
"op": "and",
"args": [
{
"op": ">=",
"args": [
{
"property": "end_datetime"
},
'2020-03-28T20:05:46+02'
]
},
{
"op": "<=",
"args": [
{
"property": "start_datetime"
},
'2020-03-28T22:06:15+02'
]
},
{
"op": "=",
"args": [
{
"property": "sar:instrument_mode"
},
'IW'
]
}
]
}
```
Optionally add spatial filtering
```python
# ext = None
ext = {'xmin': -4, 'xmax': -2, 'ymin': 6, 'ymax': 8}

if ext is not None:
arg = {
'op': 's_intersects',
'args': [
{
'property': 'geometry'
},
{
'type': 'Polygon',
'coordinates': [[[ext['xmin'], ext['ymin']],
[ext['xmin'], ext['ymax']],
[ext['xmax'], ext['ymax']],
[ext['xmax'], ext['ymin']],
[ext['xmin'], ext['ymin']]]]
}
]
}
cql2_filter['args'].append(arg)
```
Convert CQL2 filter to SQL where clause
```python
from pygeofilter.parsers.cql2_json import parse as json_parse

filter = json_parse(cql2_filter)
```

```python
from pygeofilter.backends.sql.evaluate import to_sql_where

sql_where = to_sql_where(filter, {
's1:datatake': 's1:datatake',
'datetime': 'datetime',
'sar:instrument_mode': 'sar:instrument_mode',
'end_datetime': 'end_datetime',
'start_datetime': 'start_datetime',
'geometry': 'geometry'
})
```
Create DuckDB query for a [geoparquet file](https://github.com/geopython/pygeofilter/files/14633813/20200301_20200401.zip):
Here it gets ugly because (1) the WKB column needs to be converted to `GEOMETRY` type and (2) the DuckDB-spatial implementation of `ST_GeomFromWKB` cannot read the WKB-HEX representation returned by `to_sql_where`.
```python
import re

sql_query = "SELECT * EXCLUDE geometry, ST_GeomFromWKB(geometry) AS geometry FROM '20200301_20200401.parquet' WHERE %s" % sql_where

if ext is not None:
# convert WKB blob to GEOMETRY
sql_query = sql_query.replace('ST_Intersects("geometry"', 'ST_Intersects(ST_GeomFromWKB(geometry)')

# duckdb_spatial apparently cannot yet read wkb_hex representation -> convert it back to text
spatial = ("ST_GeomFromText('POLYGON(({xmin} {ymin}, {xmin} {ymax}, "
"{xmax} {ymax}, {xmax} {ymin}, {xmin} {ymin}))')")
sql_query = re.sub(r'ST_GeomFromWKB\(x\'.*\'\)', spatial.format(**ext), sql_query)
sql_query
```
Execute the query:
```python
df = duckdb.query(sql_query)
```

I wonder, how would you do this? Do you think there is anything that could/should be modified on the `pygeofilter` end or is it entirely up to the DuckDB-spatial package? I'd appreciate any help.

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.