Geo mit Python

# Geo mit Python

---

# Über mich

### Oliver Tonnhofer

* Omniscale GmbH & Co. KG, Oldenburg
  * Open Source Entwicklung (Client/Server)
  * OpenStreetMap Kartendienste
  * MapProxy und Imposm Entwicklung und Support

---

# Was ist Python?

* Eine Programiersprache
--

* OpenSource
--

* Hat eine sehr große Community
--

* Ist einfach

---

# Einfachstes Programm

```bash
$ python hello.py
Hallo Welt!
```

---

# Interaktiv

```python
$ python
>>> print "Hallo Welt!"
Hallo Welt!
```
--

### Beispiele in den Folien
```python
>>> Befehl
Ausgabe
```

---

# iPython

```python
$ ipython
>>> import request
>>> request.<tab><tab>
[...]
requests.Timeout    requests.certs    requests.get
[...]
```
--
```python
>>> request.get?
Definition: requests.get(url, **kwargs)
Docstring:
Sends a GET request. Returns :class:`Response` object.
[...]
```
---

# Datentypen

```python
>>> l = [1000, 500, 200]
>>> l[1:]
[500, 200]
>>> l[-1]
200
```

```python
>>> d = {'jahr': 2014, 'ort': 'Berlin'}
>>> d['jahr']
2014
```

---

# Lesbar und verständlich

```python
with open(filename) as f:
    for line in f:
        print line
```

---

# *Batteries included*

### Standardbibliothek

---

# `glob`

```python
for fname in glob.glob("*.shp"):
    process_shp(fname)
```

---

# `csv`

```
id;lat;long;name;type
1;53;8;Oldenburg;city
2;52;13;Berlin;city
```

```python
with open("points.csv") as f:
    reader = csv.DictReader(f, delimiter=';')
    for r in reader:
        print r['name'], r['long'], r['lat']
```

---

# `json`

```python
with open("data.geojson") as f:
    data = json.load(f)
    for f in data['features']:
        print f['properties']
```

---

# Standardbibliothek

* noch viel mehr Pakete
* Library Reference → docs.python.org
* Python Module of the Week → pymotw.com

---

# Weitere Pakete

### Python Package Index (PyPI)

---

# Weitere Pakete

Unter Linux mit PIP:

* Robuster
* Uninstall

```shell
$ pip install neuespaket
```

Unter Windows mit Easy Install (`setuptools`):

* Binärpakete

```shell
> easy_install neuespaket
```

---

# Virtualenv

* Viele Linux Distributionen/Anwendungen verwenden Python
* Empfohlen: Pakete ausserhalb des System installieren

```python
$ virtualenv ~/meinprojekt
$ ~/meinprojekt/bin/pip install neuespaket
```

---

# Was kann ich damit anfangen?

Geodaten

* herunterladen,
* einlesen,
* transformieren,
* exportieren,
* abfragen,
* ...

---
# Requests

```python
>>> resp = requests.get('https://geo.sv.rostock.de'
...   '/download/opendata'
...   '/abfallbehaelter/abfallbehaelter.json')
```

```python
>>> if resp.ok:
...   print resp.json()
```
--

```python
{'features': [
    {'geometry': {
        'coordinates': [12.12429, 54.09293],
        'type': 'Point'},
     'properties': {
        'art': 'Abfallbehälter', 'id': 1},
     'type': 'Feature'},
```

---

# Requests

```python
  for f in resp.json()['features']:
     print f['properties']
     print f['geometry']['coordinates']
```

```json
{'art': 'Abfallbehälter', 'id': 1}
[12.124293950279831, 54.0929309062321]
{'art': 'Abfallbehälter', 'id': 2}
[12.124668944836374, 54.09288836213802]
...
```

---

# PyProj

```python
>>> p4326 = Proj(init='epsg:4326')
>>> p25832 = Proj(init='epsg:25832')
>>> transform(p4326, p25832, 12, 54)
(696620.8983966266, 5987687.710270112)
```

```python
points = []
for f in resp.json()['features']:
    long, lat = f['geometry']['coordinates']
    points.append(
        transform(p4326, p25832, long, lat)
    )
```

---

# Shapely

* Moderne API für GEOS
* Geometrien von/nach GeoJSON/WKT/WKB
* OGC *Simple Features Access Specification*
    * Buffer, Valid, Simple
    * Intersects, Crosses, Distance
    * Union, Intersection, ...

---

# Shapely

### Beispiele

```python
>>> Point(0, 0).touches(
...     LineString([(-5, -5), (5, 5)])
True
```

--
```python
>>> c1 = Point(0, 0).buffer(10)
>>> c2 = Point(10, 0).buffer(10)
>>> c1.area, c2.area
313.65484, 313.65484
>>> c1.intersection(c2).area
122.49849
>>> c1.union(c2).area
504.81120
```
---
# Shapely

### Koordinaten → Shapely Geometry

```python
point_geoms = []
for x, y in points:
    point_geoms.append(
        Point(x, y)
    )
```
--

```python
point_geoms = [Point(x, y) for x, y in points]
```

---

# Shapely

### Nach Entfernung sortieren

```python
loc = Point(702250, 5998640)
point_geoms.sort(key=lambda x: loc.distance(x))
```

```python
for p in point_geoms[:3]:
    print loc.distance(p)
```

```
31.6293983341
48.3288482691
64.938473919
```

---

# SQLAlchemy

* Datanbank API für Python
* PostgreSQL, MySQL, SQLite, etc.
* High-Level: Object-Relational-Mapper (ORM)
* Low-Level: Expression API

## GeoAlchemy

* Erweiterung für PostGIS

---

# SQLAlchemy ORM

```python
class Trash(Base):
    __tablename__ = 'trashcans'
    id = Column(Integer, primary_key=True)
    art = Column(String)
    geom = Column(Geometry('POINT', srid=25832))
```
---

# SQLAlchemy ORM

```python
engine = create_engine(
    'postgresql://localhost/gis')
session = Session(bind=engine)

for f in r.json()['features']:
    geom = WKTElement(
        shape(f['geometry']).wkt, 4326)
    t = Trash(
        art=f['properties']['art'],
        geom=func.st_transform(geom, 25832),
    )
    session.add(t)

session.commit()
```

---

---

# SQLAlchemy ORM

```python
loc = from_shape(Point(701163 6007829), 25832)

query = session.query(Trash).filter(
    Trash.geom.intersects(
        func.ST_Buffer(loc, 100)
    )
)
```
--

```python
for t in query:
    print t.art, to_shape(t.geom)
```
```
Abfallbehälter POINT (701178.013 6007874.937)
Abfallbehälter POINT (701186.835 6007849.791)
Abfallbehälter POINT (701216.373 6007762.855)
...
```

---

## SQLAlchemy Expression API

```python
poi_types = expression.case([
    (pois.c.tourism == 'hostel', 'herberge'),
    (or_(
        pois.c.amenity == 'hotel',
        pois.c.building == 'hotel',
        pois.c.tourism == 'hotel',
    ), 'hotel'),
], else_=None)

subq = select([
    pois.c.name,
    label('geometry', st_pointonsurface(pois.c.geometry)),
    label('type', poi_types])
subq = subq.where(bbox_intersects(pois.c.geometry, bbox, srs, db_srs))
subq = subq.where(pois.c.name != '')
subq = expression.alias(subq, 'pois')

q = select([
    subq.c.name, subq.c.type,
    label('geometry',
        st_astext(st_transform(st_collect(subq.c.geometry), target_srs))
        ),
])
q = q.where(subq.c.name != '')
q = q.where(subq.c.type.in_(poi_types))
q = q.group_by(subq.c.type, subq.c.name)
q = q.order_by(subq.c.type)
```

---

# Fiona

* Moderne API für OGR
* Shapefile, Mapinfo TAB, etc.
* Objekte als *"GeoJSON"*

---

# Fiona

```python
schema = {
    'geometry': 'Point',
    'properties': {'art': 'str'},
}
```

```python
with fiona.open(
    '/tmp/trash.shp', 'w',
    driver='ESRI Shapefile',
    schema=schema,
    crs=from_epsg(4326),
    ) as f:
```

```python
    f.write({
        'geometry': {
            'type': 'Point',
            'coordinates': [12, 54]
        },
        'properties': {'art': u'Abfallbehälter'},
    })
```
---

# Was gibt es noch?

---

# RTree

* API für `libspatialindex`
* Räumlicher Index
* *In welchem Polygon liegt Punkt X?*

---

# rasterio

* Moderne API für GDAL
* noch im Entwicklungsstadium
* GeoTIFF, JPEG, ECW, etc.
* Raster als NumPy

---

# NumPy

* *"the fundamental package for scientific computing"*
* Alternative zu MATLAB
* Sehr große Community mit hunderten Erweiterungen
  * SciPy
  * Matplotlib
* Algorithmen für Vektor und Raster
  * Voronoi-Diagramme, KNN, Clustering, etc. pp.

---

# Anwendungen

* QGIS: Eigene Plugins und Skripte in Python

---

---

# Anwendungen

* QGIS: Eigene Plugins und Skripte in Python
* PyWPS: WPS Prozesse in Python
* Imposm: OSM nach PostGIS

---

# Überblick

* Python: JSON, CSV, Plain-Text
--

* Requests: HTTP
--

* PyProj: Koordinaten-Transformationen
--

* Shapely: Geometrien
--

* SQLAlchemy: SQL
--

* GeoAlchemy: PostGIS
--

* Fiona: Vektordaten
--

* rasterio: Rasterdaten
--

* NumPy: Algorithmen

---

# Fazit

Python...

* bietet Pakete für jeden Zweck,
--

* ermöglicht schnelle Prototypen/Skripte,
--

* aber auch produktive Anwendung,
--

* ermöglicht lesbaren Code.

---

# Fin

* Oliver Tonnhofer
* <img src="imgs/omniscale_logo_4c_cymk.png" alt="Omniscale" style="display:inline; height: 1.4em;"/>
* `tonnhofer@omniscale.de`
* `@oltonn`