Skip to content

Commit

Permalink
chore: update after review
Browse files Browse the repository at this point in the history
Co-authored-by: Phillip Cloud <[email protected]>
  • Loading branch information
ncclementi and cpcloud committed Dec 11, 2023
1 parent 491da49 commit 9e10d77
Showing 1 changed file with 24 additions and 24 deletions.
48 changes: 24 additions & 24 deletions docs/posts/ibis-duckdb-geospatial/index.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,9 @@ categories:
- geospatial
---

Ibis now has support for some DuckDB geospatial functions (check the list [here](https://gist.github.com/ncclementi/fbc5564af709e2d7f8882821e3a8649f)).
Ibis now has support for [DuckDB geospatial functions](https://gist.github.com/ncclementi/fbc5564af709e2d7f8882821e3a8649f)!

This blogpost showcases some examples of usage of the geospatial API for the DuckDB backend. The material is inspired in
This blogpost showcases some examples of the geospatial API for the DuckDB backend. The material is inspired in
the ["DuckDB: Spatial Relationships"](https://geog-414.gishub.org/book/duckdb/07_spatial_relationships.html) lesson from
[Dr. Qiusheng Wu](https://geog-414.gishub.org/book/preface/instructor.html)'s course "Spatial Data Management" from the
Department of Geography & Sustainability at the University of Tennessee, Knoxville.
Expand Down Expand Up @@ -61,12 +61,12 @@ spatial relations.
We can start by taking a peak at the `nyc_subway_stations` table.

```{python}
subway_stations = con.table('nyc_subway_stations')
subway_stations = con.table("nyc_subway_stations")
subway_stations
```

Notice that the last column has a `geometry` type, and in this case it contains points that represent the location of
each subway station. Let's grab the entry for `'Broad St'`.
each subway station. Let's grab the entry for the Broad St subway station.

```{python}
broad_station = subway_stations.filter(subway_stations.NAME == "Broad St")
Expand Down Expand Up @@ -94,7 +94,7 @@ subway_stations.filter(_.geom.geo_equals(_.filter(_.NAME == "Broad St").geom))

### `intersect` (ST_Intersect)

Let's find in which neighborhood is "Broad Street" subway station and determine its neighborhood using the the
Let's locate the neighborhood of the "Broad Street" subway station using the
geospatial `intersect` function. The `intersect` function returns `True` if two geometries have any points in common.

```{python}
Expand All @@ -108,7 +108,7 @@ boroughs.filter(boroughs.geom.intersects(broad_station.select(broad_station.geom

### `d_within` (ST_DWithin)

We can also find what are the streets nearby (within 10 meters of) of the "Broad St" subway station using the `d_within`
We can also find the streets near (say, within 10 meters) the Broad St subway station using the `d_within`
function. The `d_within` function returns True if the geometries are within a given distance.

```{python}
Expand All @@ -128,9 +128,7 @@ In the previous query, `streets` and `broad_station` are different tables. We us
scalar subquery from a table with a single column (whose shape is scalar).
:::

To visualize the findings, we will convert the tables to Geopandas dataframes. We need to do this because,
at the moment, DuckDB geometries are not capable to store SRID information. By converting to
Geopandas we can assign a `crs` and make some nice plots.
To visualize the findings, we will convert the tables to Geopandas DataFrames.

```{python}
broad_station_gdf = broad_station.to_pandas()
Expand All @@ -144,24 +142,26 @@ streets_gdf.crs = "EPSG:26918"
```

```{python}
#visualize multiple layers
import leafmap.deckgl as leafmap
import leafmap.deckgl as leafmap # <1>
```

1. `leafmap.deckgl` allows us to visualize multiple layers

```{python}
m = leafmap.Map()
m.add_vector(broad_station_gdf, get_fill_color='blue')
m.add_vector(sts_near_broad_gdf,get_color='red', opacity=0.5)
m.add_vector(streets_gdf, get_color='grey', zoom_to_layer=False, opacity=0.3)
m.add_vector(broad_station_gdf, get_fill_color="blue")
m.add_vector(sts_near_broad_gdf, get_color="red", opacity=0.5)
m.add_vector(streets_gdf, get_color="grey", zoom_to_layer=False, opacity=0.3)
m
```

You can zoom in and out, and hover over the map to check on the streets names.

### `buffer` (ST_Buffer)

We have a table in our database with information about homicides:
Next, we'll take a look at the homicides table and showcase some
additional functionality related to polygon handling.

```{python}
homicides = con.table("nyc_homicides")
Expand All @@ -183,17 +183,17 @@ We can check the area using the `area` (`ST_Area`) function, and see that is $~
broad_station.geom.buffer(200).area()
```

To find if there were any homicides that happened in that area, we can find where the polygon resulting from adding the
To find if there were any homicides in that area, we can find where the polygon resulting from adding the
200 meters buffer to our "Broad St" station point intersects with the geometry column in our homicides table.

```{python}
h_near_broad = homicides.filter(_.geom.intersects(broad_station.select(_.geom.buffer(200)).to_array()))
h_near_broad
```

It looks like there was one homicide that happened in a radius of 200 meters from the "Broad St" station, but from this
data we can not tell in near which street happened. However, we can check if the homicide point is within a small
distance of a street and get that information.
It looks like there was one homicide within 200 meters from the "Broad St" station, but from this
data we can't tell the street near which it happened. However, we can check if the homicide point is within a small
distance of a street.

```{python}
h_street = streets.filter(_.geom.d_within(h_near_broad.select(_.geom).to_array(), 2))
Expand All @@ -215,11 +215,11 @@ h_street_gdf.crs = "EPSG:26918"
mh = leafmap.Map()
mh.add_vector(broad_station_gdf, get_fill_color='orange')
mh.add_vector(broad_station_zone, get_fill_color='orange', opacity=0.1)
mh.add_vector(h_near_broad_gdf, get_fill_color='red', opacity=0.5)
mh.add_vector(h_street_gdf, get_color='blue', opacity=0.3)
mh.add_vector(streets_gdf, get_color='grey', zoom_to_layer=False, opacity=0.2)
mh.add_vector(broad_station_gdf, get_fill_color="orange")
mh.add_vector(broad_station_zone, get_fill_color="orange", opacity=0.1)
mh.add_vector(h_near_broad_gdf, get_fill_color="red", opacity=0.5)
mh.add_vector(h_street_gdf, get_color="blue", opacity=0.3)
mh.add_vector(streets_gdf, get_color="grey", zoom_to_layer=False, opacity=0.2)
mh
```
Expand Down

0 comments on commit 9e10d77

Please sign in to comment.