GeoPandas - easy, fast and scalable geospatial analysis in Python

class: center, middle

# GeoPandas

##  Easy, fast and scalable geospatial analysis in Python

Joris Van den Bossche, PyData Paris, January 21, 2019

https://github.com/jorisvandenbossche/talks/

[@jorisvdbossche](https://twitter.com/jorisvdbossche)

---
class: center, middle
count: false

# GeoPandas

##  Extending Pandas and Dask to enable easy, fast and scalable geospatial analysis in Python

Joris Van den Bossche, PyData Paris, January 21, 2019

https://github.com/jorisvandenbossche/talks/

[@jorisvdbossche](https://twitter.com/jorisvdbossche)

---
# About me

Joris Van den Bossche

- Background: PhD bio-science engineer, air quality research
- Open source enthusiast: pandas core dev, geopandas maintainer, scikit-learn contributor
- Currently working at the Université Paris-Saclay Center for Data Science (Inria)

https://github.com/jorisvandenbossche   Twitter: [@jorisvdbossche](https://twitter.com/jorisvdbossche)

.affiliations[
  ![:scale 65%](img/logoUPSayPlusCDS_990.png)
  ![:scale 25%](img/inria-logo.png)
]

---
class: center, middle

#  Geospatial data

---

background-image: url(img/us_census_density.jpg)

---
class: center

![:scale 100%](img/inbo_bird_tracking.png)

.credits[
Bird tracking. From [slides](https://speakerdeck.com/peterdesmet/fly-my-pretties-fly-exploring-visualizing-and-publishing-bird-tracking-data) by Peter Desmet (INBO)
]

---
class: center

![:scale 90%](img/pickups_dropoffs_narrow.jpg)

.credits[
NYC taxi data (pick-up and drop-off locations), from [blog post](https://medium.com/towards-data-science/if-taxi-trips-were-fireflies-1-3-billion-nyc-taxi-trips-plotted-b34e89f96cfa) by Ravi Shekhar]

---
# Raster vs vector data

![:scale 49%](img/raster_example.png)
![:scale 49%](img/vector_example.png)

--
count: false

.right[
### -> in this talk: focus on vector data
]

.right[
### -> simple features (points, linestrings, polygons) with attributes
]

???

Two major families of geospatial data

raster: grid based (topology lacking, difficult to link to tabular data)
vector: coordinate based objects, topological

here: vector

vector -> common abstraction model in many software
Open Geospatial consortium standard

---

# OGC Simple Features

Point, LineString, Polygon

.center[
![:scale 60%](img/simple_features.svg)
]

... and MultiPoint, MultiLineString, MultiPolygon, GeometryCollection

--
count: false

Attributes : each vector feature can have a record in attribute table

???

and that is where geopandas comes into play

but before talking about geopandas, first a bit more general about open source geospatial software

# geospatial software

This presentation: in python

but everything I will present -> builds upon widely used open source libraries

---
class: middle, center

# Open source geospatial software

.center[
![:scale 70%](img/Open_Source_Geospatial_Foundation.svg)
]

???

Open Source Geospatial Foundation

OSGeo was created to support the collaborative development of open source geospatial software, and promote its widespread use.

---

# GDAL / OGR

### Geospatial Data Abstraction Library.

* The swiss army knife for geospatial.
* Read and write Raster (GDAL) and Vector (OGR) datasets
* More than 200 (mainly) geospatial formats and protocols.

.center[
![:scale 100%](img/gdal_formats)
]

.credits[
Slide from "GDAL 2.2 What's new?" by Even Rouault (CC BY-SA)
]

???

GDAL is a translator library for raster and vector geospatial data formats. As a library, it presents a single raster abstract data model and single vector abstract data model to the calling application for all supported formats. It also comes with a variety of useful command line utilities for data translation and processing.

<!--
# GDAL / OGR

### Widely used (FOSS & proprietary)

.center[
![:scale 100%](img/gdal_users)
]

.credits[
Slide from "GDAL 2.2 What's new?" by Even Rouault (CC BY-SA)
] -->

---

# GEOS

## Geometry Engine Open Source

* C/C++ port of a subset of Java Topology Suite (JTS)
* Most widely used geospatial C++ geometry library
* Implements geometry objects (simple features), spatial predicate functions and spatial operations

Used under the hood by many applications (QGIS, PostGIS, MapServer, GRASS, GeoDjango, ...)

[geos.osgeo.org](http://geos.osgeo.org)

<!-- ---
# PROJ.4

C library for performing conversions between cartographic projections.

[http://proj4.org/](http://proj4.org/) -->

---

# Python geospatial packages

--
count:false

Interfaces to widely used libraries:

- Python bindings to GDAL/OGR (`from osgeo import gdal, ogr`)
- [`pyproj`](https://jswhit.github.io/pyproj/): python interface to PROJ.4.
- Pythonic binding to GDAL/OGR:
  - [`rasterio`](https://mapbox.github.io/rasterio/) for GDAL
  - [`fiona`](http://toblerity.org/fiona/README.html) for OGR
- [`shapely`](https://shapely.readthedocs.io/en/latest/): python package based on GEOS.

---

# Shapely

Python package for the manipulation and analysis of geometric objects

Pythonic interface to GEOS

--
count:false

.mmedium[
```python
>>> from shapely.geometry import Point, LineString, Polygon

>>> point = Point(1, 1)
>>> line = LineString([(0, 0), (1, 2), (2, 2)])
>>> poly = line.buffer(1)
```
]

--
count: false

Nice interface to GEOS, but: single objects, no attributes

???

# Shapely

typical predicates and operations

(images from shapely docs)

---

![:scale 100%](img/pandas_logo.svg)

One of the packages driving the growing popularity of Python for data science, machine learning and academic research

* High-performance, easy-to-use data structures and tools
* Suited for tabular data (e.g. columnar data, spread-sheets, database tables)

```python
import pandas as pd

df = pd.read_csv("myfile.csv")

subset = df[df['value'] > 0]
subset.groupby('key').mean()
```

---
class: middle, center

# GeoPandas

# .darkred[.fat[Easy]], fast and scalable geospatial analysis in Python

---

# GeoPandas

Make working with geospatial data in python easier

* Started by Kelsey Jordahl in 2013
* Extends the pandas data analysis library to work with geographic objects and spatial operations
* Combines the power of whole ecosystem of (geo) tools (pandas, geos, shapely, gdal, fiona, pyproj, rtree, ...)
* Bridge between geospatial packages and the scientific / data science stack

Documentation: http://geopandas.readthedocs.io/

???

make working with geospatial data like working with any other kind of data in python
(data stack, numpy, pandas and other tools around those)
analysis for which you otherwise would need desktop GIS applications (QGIS, ArcGIS) or geospatial databases (PostGIS)

makes pandas objects geometry aware

---
class: middle, center

# Demo time!

See [static version](https://nbviewer.jupyter.org/github/jorisvandenbossche/geopandas-demo/blob/master/geopandas_demo.ipynb)

---

# Summary

* Read and write variety of formats (fiona, GDAL/OGR)
* Familiar manipulation of the attributes (pandas dataframe)
* Element-wise spatial predicates (intersects, within, ...) and operations (intersection, union, difference, ..) (shapely)
* Re-project your data (pyproj)
* Quickly visualize the geometries (matplotlib, descartes)
* More advanced spatial operations: spatial joins and overlays (rtree)

--
count:false

**-> Interactive exploration and analysis of geospatial data**

---

# Ecosystem

[geoplot](http://www.residentmar.io/geoplot/index.html) (high-level geospatial visualization), [cartopy](http://scitools.org.uk/cartopy/) (projection aware cartographic library)

[folium](https://github.com/python-visualization/folium), [ipyleaflet](https://github.com/jupyter-widgets/ipyleaflet) (Leaflet.js maps)

[OSMnx](http://geoffboeing.com/2016/11/osmnx-python-street-networks/) (python for street networks)

[PySAL](http://pysal.readthedocs.io/en/latest/index.html) (Python Spatial Analysis Library)

[rasterio](https://mapbox.github.io/rasterio/) (working with geospatial raster data)

...

---
class: middle, center

# Intermezzo
# --
# Extending Pandas

---
# Extending Pandas

Composition versus inheritance

--
count: false

Subclassing pandas.DataFrame

---
# Subclassing DataFrame

```python
class GeoDataFrame(pd.DataFrame):

_metadata = ['crs', '_geometry_column_name']

def __init__(self, ...): ...

@property
    def _constructor(self):
        return GeoDataFrame

# override some methods

def plot(self, ...):
        ...

# add additional geo-specific properties and methods

@property
    def geometry(self):
        ...

```

???

Subclassing has some disadvantages:

* overriding methods can be confusing / surprising
* can loose your type (conversion to basic DataFrame, you loose functionality)

---
# Extending Pandas

Composition versus inheritance

Subclassing pandas.DataFrame

New methods since pandas 0.23:

* Registering custom accessors (for DataFrame, Series and Index)
* Extension Arrays to add custom data type support to pandas

See more in the [Extending Pandas docs](http://pandas.pydata.org/pandas-docs/stable/extending.html)

---

# Extending pandas data types

A new Extension Array interface: define your own array-like and tell pandas how to work with it.

--
count: false

Examples with [cyberpandas](https://cyberpandas.readthedocs.io/en/latest/index.html) and [GeoPandas](http://geopandas.readthedocs.io/en/latest/):

```python
>>> df
   ID    addresses           location
0   0  192.168.1.1  POINT (48.8  2.3)
1   1      0.0.0.0   POINT (51.2 4.4)

>>> df.dtypes
ID              int64
addresses          ip
location     geometry
dtype: object
```

---
count: false

# Extending pandas data types

A new Extension Array interface: define your own array-like and tell pandas how to work with it.

Examples with [cyberpandas](https://cyberpandas.readthedocs.io/en/latest/index.html) and [GeoPandas](http://geopandas.readthedocs.io/en/latest/):

```python
>>> df
   ID    addresses           location
0   0  192.168.1.1  POINT (48.8  2.3)
1   1      0.0.0.0   POINT (51.2 4.4)

>>> df.dtypes
ID              int64
*addresses          ip
*location     geometry
dtype: object
```

Each of those libraries provide additional type-specific functionality.

---

# Extending pandas data types

Motivation

* Pandas historically bound to NumPy type system and its limitations (missing data non-float dtypes, non-numeric dtypes as categorical, datetime with timezone, variable-length strings, ...)
* The custom types in pandas (Categorical, Period, Interval, ...) caused complex internals

---

# Extending pandas data types

`ExtensionDtype`

* Name and what type of scalars?

`ExtensionArray`

* Class which does the actual "heavy lifting" (storage, basic array ops)
* Some required methods/attributes that pandas needs
* Free to have additional functionality, no restriction on construction
* Limited to one dimension, though may be backed by 0..n arrays

A Series is a container for an “array-like” thing

---
class: middle, center

# Demo time!

See [static version](https://nbviewer.jupyter.org/github/jorisvandenbossche/talks/blob/master/2019_PyDataParis_geopandas_extending_pandas/pandas-extension-array-demo.ipynb)

---
class: middle, center

# GeoPandas

#  Easy, .darkred[.fat[fast]] and scalable geospatial analysis in Python

---

# However ...

---
count: false

# However ... it can be slow

Timings for basic `within` and `distance` operation on 100 000 points:

.medium[
```python
s.within(polygon)
s.distance(polygon)
```
]

![:scale 49%](img/timings_within.png)
![:scale 49%](img/timings_distance.png)

---

# Comparison with PostGIS

![:scale 49%](img/timings_distance2.png)
![:scale 49%](img/timings_sjoin.png)

Disclaimer: dummy benchmark, and I am not a PostGIS expert!

Example from [Boundless tutorial](http://workshops.boundlessgeo.com/postgis-intro/) (CC BY SA)

---

# Comparison with PostGIS

.small[
```sql
-- What is the population and racial make-up of the neighborhoods of Manhattan?
SELECT
  neighborhoods.name AS neighborhood_name, Sum(census.popn_total) AS population,
  100.0 * Sum(census.popn_white) / NULLIF(Sum(census.popn_total),0) AS white_pct,
  100.0 * Sum(census.popn_black) / NULLIF(Sum(census.popn_total),0) AS black_pct
FROM nyc_neighborhoods AS neighborhoods
JOIN nyc_census_blocks AS census
ON ST_Intersects(neighborhoods.geom, census.geom)
GROUP BY neighborhoods.name
ORDER BY white_pct DESC;
```
]

.small[
```python
res = geopandas.sjoin(nyc_neighborhoods, nyc_census_blocks, op='intersects')
res = res.groupby('NAME')[['POPN_TOTAL', 'POPN_WHITE', 'POPN_BLACK']].sum()
res['POPN_BLACK'] = res['POPN_BLACK'] / res['POPN_TOTAL'] * 100
res['POPN_WHITE'] = res['POPN_WHITE'] / res['POPN_TOTAL'] * 100
res.sort_values('POPN_WHITE', ascending=False)
```
]

Disclaimer: dummy benchmark, and I am not a PostGIS expert!

Example from [Boundless tutorial](http://workshops.boundlessgeo.com/postgis-intro/) (CC BY SA)

---

# Why is GeoPandas slower?

- GeoPandas stores custom Python objects in arrays
- For operations, it iterates through those objects
- Those Python objects each call the GEOS C operation

.center[
![:scale 60%](img/geopandas-shapely-1.svg)
]

---

# New version in development

.center[
![:scale 70%](img/geopandas-shapely-2.svg)
]

---

# New version in development

Old object-based python loops:

```python
def distance(self, other):
    result = [geom.distance(other) for geom in self.geometry]
    return pd.Series(result)

```

---

# New version in development

New C-level loop on GEOS objects:

```python
cpdef distance(self, other):
    cdef int n = self.size
    cdef double[:] out = np.empty(n, dtype=np.float64)
    cdef GEOSGeometry *left_geom
    cdef GEOSGeometry *right_geom = other.__geom__
    geometries = self._geometry_array

with nogil:
        for idx in xrange(n):
            left_geom = <GEOSGeometry *> geometries[idx]
            if left_geom != NULL:
                out[idx] = GEOSDistance_r(
                    left_geom, some_point.__geom)
            else:
                out[idx] = NaN
    return out
```

---

# New version in development

.center[
![:scale 70%](img/geopandas-shapely-2.svg)
]

Remove python overhead by only storing pointers to C GEOS objects and iterating in C

TL;DR: same API, but better performance and less memory use

Many thanks to Matthew Rocklin (Anaconda, Inc.) for his work!

---

# New timings

![:scale 40%](img/timings_within.png)
![:scale 40%](img/timings_distance.png)

![:scale 40%](img/timings_distance2.png)
![:scale 40%](img/timings_sjoin.png)

---
count:false

# New timings

![:scale 40%](img/timings_within_all.png)
![:scale 40%](img/timings_distance_all.png)

![:scale 40%](img/timings_distance2_all.png)
![:scale 40%](img/timings_sjoin_all.png)

---
# Sounds interesting?

Blogpost of me and Matthew with more background:

* http://matthewrocklin.com/blog/work/2017/09/21/accelerating-geopandas-1
* https://jorisvandenbossche.github.io/blog/2017/09/19/geopandas-cython/

Try out development version (binary builds):

.medium[
```bash
conda install --channel conda-forge/label/dev geopandas
```
]

.small[
```bash
conda create -n geo_test python=3.6 geopandas --channel conda-forge/label/dev
```
]

---
class: middle, center

# GeoPandas

#  Easy, fast and .darkred[.fat[scalable]] geospatial analysis in Python

???

Python has a fast and pragmatic data science ecosystem
... restricted to in-memory and a single core

---
.center[
![:scale 55%](img/dask_horizontal.svg)
]
## A flexible library for parallelism

--
count: false

* A parallel computing framework, written in pure Python
* Lets you work on larger-than-memory datasets
* That leverages the excellent Python ecosystem
* Using blocked algorithms and task scheduling

http://dask.pydata.org/

---

# An experiment with taxi data

[Ravi Shekhar](http://people.earth.yale.edu/profile/ravi-shekhar/about) published a blogpost [Geospatial Operations at Scale with Dask and GeoPandas](https://medium.com/towards-data-science/geospatial-operations-at-scale-with-dask-and-geopandas-4d92d00eb7e8) in which he counted the number of rides originating from each of the official taxi zones of New York City

Matthew Rocklin re-ran the experiment with the in-development version: 3h -> 8min ([see his blogpost](http://matthewrocklin.com/blog/work/2017/09/21/accelerating-geopandas-1))

[dask-geopandas](https://github.com/mrocklin/dask-geopandas): experimental library with parallelized geospatial operations and joins

--
count: false

## Demo time!

See [static version](http://nbviewer.jupyter.org/gist/jorisvandenbossche/67be41a246c1281d7046b31690988321)

--
count: false

Experimental, but ready to experiment with!

<!-- # Want to try this out?

Feedback and polishing needed!

See https://github.com/geopandas/geopandas/issues/473

Build new version (needs GEOS and shapely already installed):

* clone https://github.com/geopandas/geopandas
* checkout 'geopandas-cython' branch
* `make inplace`

# GeoPandas - dask bridge

https://github.com/mrocklin/dask-geopandas

### **Experimental** library with parallelized geospatial operations and joins

Many areas for improvement:

- higher coverage of functionality (overlay, ...)
- better use of spatial indexing (e.g. persist data on disk partitioned by region for more efficient intersects, joins, ..)
- better integration with parquet
- ...

Experimental, but ready to be tried out!
 -->

---

# Contributions welcome!

### Feedback, bug reports, feature requests

### Extend functionality

### Contribute code: https://github.com/geopandas/geopandas

---
class: middle

http://geopandas.readthedocs.io

# Thanks for listening!

## Thanks to all contributors!

## Those slides:

- https://github.com/jorisvandenbossche/talks/

http://geopandas.readthedocs.io