Spatial Data Analysis with DuckDB
Tags: python, duckdb, spatial, opendata
Published: 2023-08-06 19:18:20

## Introduction
In this article, I will explain how to access point of interest (POI) data published as open data by the Overture Maps Foundation with DuckDB and how to perform spatial analysis.

## What is DuckDB?
DuckDB is designed to support analytical query workloads, also known as Online analytical processing (OLAP). It includes a columnar-vectorized query execution engine. This is more performant than traditional systems such as PostgreSQL, MySQL, or SQLite, which process each row sequentially. There are many plugins available. You can easily transfer your data in environments such as Amazon S3, Google Cloud Storage, postgresql using plugins. You can perform spatial analysis by installing the Spatial plugin.

**Environment Setup**

```
pip install duckdb==0.8.1
```
**Creating a database**

```python
import duckdb
db = duckdb.connect("data.db")
```

**Installing plug-ins for data access and spatial analysis**
We install the "spatial" plugin to perform spatial analysis.
We are installing the "httpfs" plugin to access POI data in Amazon S3. Then we define the region as "us-west-2".

```python
db.sql("""
INSTALL spatial;
INSTALL httpfs;
LOAD spatial;
LOAD httpfs;
SET s3_region='us-west-2';
""")
```

We create a table in the database and transfer the data in parquet format. 59 million rows of data were transferred to the database in about 34 minutes.

```python
db.sql("""
  create table places as 
  select * from read_parquet('s3://overturemaps-us-west-2/release/2023-07-26-alpha.0/theme=places/type=*/*')
""")
```

```python
db.sql("""
select count(*) as count from places
""").show()
```

![count places](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/5v8spvgyzekldsc9i5pq.JPG)

```python
db.sql("""
    select * from places limit 5
""").show()
```

![limit 5](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/zywl6kbqxsp7s9fkjiao.JPG)

You can find the diagram of the POI data [here](https://medium.com/r/?url=https%3A%2F%2Fdocs.overturemaps.org%2Freference%2Fplaces%2Fplace%2F). There are columns in the data that we need to preprocess.

```python
db.sql("""
    select names, categories, confidence,brand,addresses from places limit 5
""").show()
```

![columns](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ihs1i7drvcpwopqjp08n.JPG)

For example, in order to find out which category it is in the categories column, we need to get the information from the data held in the "struct" type. You can review the document to learn about DuckDB data types.
For example, to extract which country you are located in the "Addresses" column:

```python
db.sql("""
    select replace(json_extract(CAST(addresses AS JSON), '$[0].country')::varchar,'"','') as country from places limit 5
""").show()
```

![country names](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/2kxa9p22szkg43zll46t.JPG)
After creating a column called “country” to extract country short names, we add the extracted data.

```python
db.sql("""ALTER TABLE places ADD COLUMN country VARCHAR;
       update places set country = replace(json_extract(CAST(places.addresses AS JSON), '$[0].country')::varchar,'"','')
    
""")
```
We run the following query to add the POI data in Turkey to a separate table and obtain the address, category, name, geometry information.

```python
db.sql("""
       create or replace table turkey_places as (
              select
                     replace(json_extract(places.addresses::json,'$[0].locality'),'"','')::varchar as locality,
                     replace(json_extract(places.addresses::json,'$[0].region'),'"','')::varchar as region,
                     replace(json_extract(places.addresses::json,'$[0].postcode'),'"','')::varchar as postcode,
                     replace(json_extract(places.addresses::json,'$[0].freeform'),'"','')::varchar as freeform,

                     categories.main as categories_main,

                     replace(json_extract(places.names::json,'$.common[0].value'),'"','')::varchar as names,
                     confidence,
                     bbox,
                     st_transform(st_point(st_y(st_geomfromwkb(geometry)),st_x(st_geomfromwkb(geometry))),'EPSG:4326','EPSG:3857') as geom
                     
       
              from places 
                     where country ='TR' 
       )
       
    
""")
```

Created table:

```python
db.sql("""
select * from turkey_places limit 5
""").df()
```

![turkey places](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ywbu433tmzjmekxpv4a6.JPG)

As an example, I will examine the POI data in Istanbul. I created two tables to obtain the POI points located within 500 m of the designated park points.

```python
db.sql("""
    create or replace table park_ist as (
        select * from turkey_places where locality = 'İstanbul' and categories_main='park'   
    );

    create or replace table poi_ist as (
        select * from turkey_places where locality = 'İstanbul' and categories_main <> 'park'
    )
       
""")
```
Number of POIs in Istanbul:

```python
db.sql(
    """
select count(*) from poi_ist

"""
)
'''
count
  181959
'''
```

Number of POIs designated as Parks in Istanbul:

```python
db.sql(
    """
select count(*) from park_ist

""")

'''
count
  492
'''
```

To query the POI points within 500 m of the points included in the park category:

```python
df = db.sql("""
              select  poi_ist.region as poi_ist_region,poi_ist.freeform as poi_ist_freeform,poi_ist.categories_main as poi_ist_categori ,
              park_ist.categories_main as park_categori , park_ist.names as park_names, park_ist.freeform as park_ist_freeform,

              st_distance(poi_ist.geom,park_ist.geom) as dist,
              ST_AsText(poi_ist.geom) as geom,
              ST_AsText(park_ist.geom) as geom2

              from poi_ist, park_ist 

              where ST_DWithin(poi_ist.geom, park_ist.geom,500) 
              

              
              
       
       """).to_df()

gdf = gpd.GeoDataFrame(df,geometry= gpd.GeoSeries.from_wkt(df['geom']),crs="EPSG:3857")
gdf.to_file("export/poi.geojson",driver="GeoJSON")
```

I converted the resulting table to pandas dataframe format. I then saved the table as geojson format using geopandas. I visualized it by category info with QGIS.

![visualize](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/jx1jv2zosvojdwbmrjp2.JPG)

POI data can be examined in more detail. I wanted to tell you how to work with spatial data with DuckDB. I hope it was useful. Hope to see you in the next article.




**Source:**
[Overture Maps Foundation](https://docs.overturemaps.org)
[POI Data](https://docs.overturemaps.org/reference/places/place/)
[DuckDB](https://duckdb.org/)
[DuckDB Spatial extension](https://duckdb.org/docs/extensions/spatial)
[DuckDB HTTPFS extension](https://duckdb.org/docs/extensions/httpfs)