Introduction to GIS scripting

class: center, middle

# Introduction to GIS scripting
## Python for GIS

INBO workshop, 9/11 May 2017

Stijn Van Hoey

https://github.com/stijnvanhoey/course_gis_scripting

---
class: center, middle

# Recap

---
class: left, middle

* `Packages` provide new objects and methods

```python
import datetime
import numpy as np
import pandas as pd
```

* `methods/functions` act on `objects` to do something -> `()`

```python
np.sum(my_array)
my_df.transpose()
```

* These objects have `characteristics/attributes`

```python
my_array.size
my_df.columns
```

* Some of these objects can be `sliced` -> `[]`

```python
my_array[1:10:2]
my_df["column_name"]
```

---

# Afternoon session

* An introduction example
* The GIS-ecosystem in Python
* Vectors: the essentials
* Raster suggestions
* Desktop GIS interaction

---
class: center, middle

# Introduction showcase

`03-geopandas-appetizer.ipynb`

---
class: center, middle

# Python's GIS ecosystem

#### ## With a subjective flavour/selection
 
---
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-1.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-2.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-3.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-4.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-5.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-14.png)
background-size: cover

---
class: center, middle

![:scale 10%](./img/gdal_logo.png)

[GDAL (Geospatial Data Abstraction Library)](http://www.gdal.org/) is the 
open source Swiss Army knife for raster and vector geospatial data formats.

Traditionally `GDAL` used to design the raster part of the library, and `OGR` the vector part for Simple Features. Starting with GDAL 2.0, both sides have been more integrated.

---

GDAL's python bindings expose most of the functionality of GDAL.

```python
from osgeo import gdal
from osgeo import ogr
```

### But...

GDAL python bindings are not very "pythonic"

```python
geo = gdal.Open(raster_file)
drv = geo.GetDriver()
print(drv.GetMetadataItem('DMD_LONGNAME'))
```

.center[
![:scale 80%](./img/monty_python.jpg)
]

---
class: left, middle

### More pythonic libraries for reading GIS data

* OGR (vector) -> Fiona
* GDAL (raster) -> Rasterio

---

## fiona

Fiona is a Python package for reading (and writing) `vector data`.

Fiona provides `Python objects` (e.g. a dictionary for each record) to geospatial data in various formats.

```python
> with fiona.open('../data/deelbekkens/Deelbekken.shp') as deelbekkens:
     feature = next(iter(deelbekkens)) # Just one checking the first
     print("Bekken: ", feature['properties']['BEKNAAM'])
     print("Vectortype: ", feature['geometry']['type'])
Bekken:  Boven-Scheldebekken
Vectortype:  Polygon
```

---

## shapely

Shapely is a Python library for geometric operations (**GEOS** library).

Shapely can perform:

* geometry validation
* geometry creation (e.g. collections)
* **geometry operations**

`difference`, `intersection`, `union`,...
.center[
![:scale 80%](./img/shapely_intersection.png)
]

---
class: middle

**Compare geometries**

* `object.almost_equals(other[, decimal=6])`
* `object.contains(other)`
* `object.crosses(other)`
* `object.disjoint(other)`
* `object.equals(other)`
* `object.intersects(other)`
* `object.touches(other)`
* `object.within(other)`

---

## pyproj
### handling Spatial Reference Systems (SRS)

Examples:
* `EPSG:4326` -> latitude, longitude in WGS-84 coordinate system
* `EPSG:31370` -> Belge 1972 / Belgian Lambert 72

A number of formats exist to represent the SRS: 
* WKT
* PROJ.4
* EPSG:xxxx

---
class: center, middle

# Vectors: the essentials

`04-gis-python-vectors.ipynb`

---
class: center, middle

## The curious case of GDAL

---

## GDAL functions: CMD

Available on the command line... with [documentation](http://gdal.org/1.11/gdalwarp.html)...

.center[
![:scale 70%](./img/ogr2ogr_cmd.png)
]

This [online cheat-sheet](https://github.com/dwtkns/gdal-cheat-sheet) is a great resource!

---

## GDAL functions: QGIS

* Add [gdaltools plugin](https://docs.qgis.org/2.2/en/docs/user_manual/plugins/plugins_gdaltools.html) to QGis
.center[
![:scale 70%](./img/gdaltools.png)
]

---

## GDAL functions: QGIS

* Build the command using the interface and copy paste your command for later usage...

.center[
![:scale 70%](./img/gdaltools_cmd.png)
]

---

## GDAL functions: Python

* Import the GDAL library and use it within a script

```python
from osgeo import gdal

#...
```

* Write a Python function to call a specific GDAL command

```python
# No 'from osgeo import gdal' (!)
subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, 
                     '-crop_to_cutline', '-overwrite'])
```

---
class: center, middle

# call GDAL

`05-the-power-of-gdal.ipynb`

---
## Rasterio

(*All credits to [Kelsey Jordahl](https://github.com/kjordahl/SciPy-Tutorial-2015) for the example*)

---

### Rasterio python API

```python
with rasterio.open('manhattan.tif') as f:
    img = f.read(1)
imshow(img, cmap='gray')
```

.center[
![:scale 80%](./img/manhattan_band_1.png)
]

---

### Rasterize vector features

Using GeoPandas to generate a vector mask and applying it to a raster image with rasterio ([example](https://github.com/kjordahl/SciPy-Tutorial-2015/blob/master/examples/rasterio_mask.py)).

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

---

### Rasterize vector features

Using GeoPandas to generate a vector mask and applying it to a raster image with rasterio ([example](https://github.com/kjordahl/SciPy-Tutorial-2015/blob/master/examples/rasterio_mask.py)).

```python
from rasterio.features import rasterize
mask = rasterize([poly], transform=src.transform, out_shape=src.shape)
```

.center[
![:scale 75%](./img/mask.png)
]

---
class: center, middle

# Raster data in Python

(*only for illustration purposes*)

`06-gis-python-rasters.ipynb`

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-6.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-7.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-8.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-9.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-10.png)
background-size: cover

---
count: false
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-11.png)
background-size: cover

---

# arcpy

```python
import arcpy
```

The `arcpy` Package provides a [set of `objects` and `functions`](http://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-functions/alphabetical-list-of-arcpy-functions.htm)

--
 
For example, the interaction with Numpy:

```python
arcpy.env.workspace = os.getcwd() + "\\demo.gdb" # connect to gdb

arr = arcpy.da.TableToNumPyArray("canton", "pop_2008")
print("Sum : {:,}".format(arr['pop_2008'].sum()))
```

--

<div class="alert alert-info">
Hence, the same as we've been doing the whole day... (except of the licen$e)
</div>

You can also [export](http://desktop.arcgis.com/en/arcmap/10.3/analyze/modelbuilder/exporting-a-model-to-a-python-script.htm) a model builder setup as a Python script and start working on it...

---
class: center, middle

# QGis scripts in Python

`split vector layer by attribute`-demo

---
class: center, middle

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

---
count: false
class: center, middle

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

---
count: false
class: center, middle

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

---
count: false
class: center, middle

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

---
class: center, middle

# QGis scripts in Python

`split vector layer by attribute`-demo

Check the information [here](https://docs.qgis.org/2.8/en/docs/user_manual/processing/scripts.html) and [here](http://docs.qgis.org/2.0/en/docs/user_manual/processing/console.html) to write your own QGIS scripts.

---
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-11.png)
background-size: cover

---
class: center, middle, bgheader
background-image: url(./img/StijnVH-ecosystem-gis-11_alpha.png)
background-size: cover

<div class="alert alert-success">
Don't be overwhelmed, enlarge your toolbox function by function...
</div>

---
class: center, middle

... and help each other!

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

---
class: center, middle

## Good luck!

![](http://www.backwaterreptiles.com/images/pythons/angolan-python-for-sale.jpg)