VectorPython.Rmd

---
author: "Jan Verbesselt, Jorge Mendes de Jesus, Aldo Bergsma, Dainius Masiliunas"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  knitrBootstrap::bootstrap_document:
    title: "Week 3: Python for geo-scripting"
    theme: "simplex"
    highlight: Tomorrow Night Bright
    menu: FALSE
    theme.chooser: TRUE
    highlight.chooser: TRUE
---

# [WUR Geoscripting](https://geoscripting-wur.github.io/) ![WUR logo](http://www.wageningenur.nl/upload/f9a0b5f2-15c5-4b84-9e14-20ef02f5e265_wur-logo.png)

# Vector data handling with Python

_Jan Verbesselt, Jorge Mendes de Jesus, Aldo Bergsma, Dainius Masiliunas_ - `r format(Sys.time(), '%d %B, %Y')`

## GDAL/OGR using Python

By now you should know what the following abbreviations stand for, and what can you do with the software:

* GDAL
* OGR

Whereas the following may be new to you: 

* OSR: spatial reference (projections!)


## OGR and the data drivers

OGR supports many different vector formats:

* ESRI formats such as shapefiles, personal geodatabases and ArcSDE
* Other software such as MapInfo, GRASS, Microstation
* Open formats such TIGER/Line, SDTS, GML, KML
* Databases such as MySQL, PostgreSQL, Oracle Spatial, etc.
* Webservices such as Web Feature Service 1.0,1.1,2.0 (WFS output as GML) 

**OGR Data drivers:** A driver is an object that knows how to interact with a certain data type (e.g. a shapefile).
You need to know the right driver to read or write data. Some drivers might only be able to read content, but the majority can read and write.

Start a terminal and see if you can load osgeo in Python.
```{r, engine='bash', eval = FALSE}
python
```

```{r, engine='python', eval = FALSE}
try:
  from osgeo import ogr
except:
  import ogr
```

From the terminal you can also type the following to obtain an overview of supported formats: 
```{r, eval=FALSE,engine='bash'}
echo "Terminal command to obtain an overview of supported formats"
ogrinfo --formats
```

[More info about the different OGR formats](http://www.gdal.org/ogr_formats.html).

## Setting up environment
Before we continue with scripting, you need to setup your Python environment. In the previous Python lesson you learned about Conda and Jupyter Notebooks. We will use both again today. Check if you have made any conda environment named geoscripting.

```{r, eval=FALSE,engine='bash'}
## Display all your conda environments
conda info --envs
```

If you don't have a conda environment called 'geoscripting', create it first.

```{r, eval=FALSE,engine='bash'}
## Create conda environment called geoscripting
conda create --name geoscripting python=3 numpy jupyter
```

Then install todays modules in your conda environment and start your Jupyter Notebook.

```{r, eval=FALSE,engine='bash'}
## Install gdal and folium
conda install --name geoscripting --channel conda-forge gdal folium
source activate geoscripting
jupyter notebook
```

Open a new Ipython Notebook in an appropriate location. Output of the Python scripts from Jupyter Notebook will be saved in the repository where the .ipynb is stored.

## Points
Now we continue with the tutorial. Try the script below and find out what "WKT" means. 

```{r, engine='python'}
from osgeo import ogr
# Create a point geometry
wkt = "POINT (173914.00 441864.00)"
pt = ogr.CreateGeometryFromWkt(wkt)
print(pt)
# print(help(osgeo.ogr))
```

WKT (Well Known Text) is a sort of markup language that describes spatial information in a clean text format, there is an equivalent in binary format (easier for computers to process and more efficient for data transfer).

Explanation of WKT: [Wikipedia](https://en.wikipedia.org/wiki/Well-known_text)


## Defining a spatial reference

A spatial reference defines: an ellipsoid, a datum using that ellipsoid, and either a geocentric, geographic or projected coordinate system.  

Therefore spatial references come in multiple forms according to our data needs. There is even a search engine for it: [http://spatialreference.org](http://spatialreference.org/ref/epsg/wgs-84/).

The following example is based on WGS84, see: [http://spatialreference.org/ref/epsg/wgs-84/](http://spatialreference.org/ref/epsg/wgs-84/)

```{r, engine='python'}
from osgeo import osr
## Spatial reference
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(4326)  # from EPSG - Lat/long
```

## Reproject a point

```{r, engine = 'python', eval = FALSE}
## lat/long definition
source = osr.SpatialReference()
source.ImportFromEPSG(4326)

# http://spatialreference.org/ref/sr-org/6781/
# http://spatialreference.org/ref/epsg/28992/
target = osr.SpatialReference()
target.ImportFromEPSG(28992)

# transform from lat/long to Dutch coord. system
transform = osr.CoordinateTransformation(source, target)
point = ogr.CreateGeometryFromWkt("POINT (5.6660 51.9872)") #ogr uses lon, lat
point.Transform(transform)
print(point.ExportToWkt())
```

Are you still keeping track of your data and coordinate system? 

> **Question 1**: On the previous example where WKT is printed out, in which coordinate system is the WKT represented?

## Create a shape file

**Concept:**

* A point is a type of geometry stored as a feature.
* A layer can have many features.
* A datasource can have many layers.
* The driver saves the datasource in a specific [format](http://www.gdal.org/ogr_formats.html).

This can be visualised as a pyramid of elements:

```{r, eval=FALSE}
Driver
    Datasource
        Layer
            Feature
                Geometry
                    Point
```

Now we know how a file is made, you can continue to make a file:

* Set spatial reference
* Create shape file
* Create layer
* Create point
* Put point as a geometry inside a feature
* Put feature in a layer
* Flush

```{r, engine = 'python', eval=FALSE}
## Is the ESRI Shapefile driver available?
driverName = "ESRI Shapefile"
drv = ogr.GetDriverByName(driverName)
if drv is None:
    print("%s driver not available.\n" % driverName)
else:
    print("%s driver IS available.\n" % driverName)
```

After loading modules and setting the correct driver, we can create a new data source for a shapefile. Have a look that the layer does not exist in your data folder.

```{r, engine = 'python', eval=FALSE}
## Make name for file and layer
fn = "points.shp"
layername = "anewlayer"

## Create shape file
ds = drv.CreateDataSource(fn)
print(ds.GetRefCount())
```

We created the datasource, but didn't make the file yet. Let's continue to add information to the shapefile. First add a coordinate reference system (CRS) or spatial reference.

```{r, engine = 'python', eval=FALSE}
# Set spatial reference by importing a projection
spatialReference = osr.SpatialReference()
spatialReference.ImportFromProj4('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')

# Optionally you could enter epsg code to define the crs
# spatialReference.ImportFromEPSG(4326)
```

Now we are ready to create the layer and create files.

```{r, engine = 'python', eval=FALSE}
## Create Layer
layer=ds.CreateLayer(layername, spatialReference, ogr.wkbPoint)
```

When you check your data folder, you will see that the file has been created! From now on it is not possible anymore to CreateDataSource with the same name in your work directory until you remove the name.shp, name.shx and name.dbf file.

> **Question 2**: What is the geometry type of the created layer? 

> **Question 3**: What does WKB mean? (hint: think about WKT)

```{r, engine = 'python', eval=FALSE}
## Check the extent of the layer
print(layer.GetExtent())
```

The layer does not have geometry features yet. Ok we will continue downwards on the pyramid and start building the feature and geometry.

```{r, engine = 'python', eval=FALSE}
## Create two point objects
point1 = ogr.Geometry(ogr.wkbPoint)
point2 = ogr.Geometry(ogr.wkbPoint)

## Add geometry
# SetPoint(self, int point, double x, double y, double z = 0)
point1.SetPoint(0,4.897070,52.377956) 
point2.SetPoint(0,5.104480,52.092876)

## Feature is defined from properties of the layer, e.g:
layerDefinition = layer.GetLayerDefn()
feature1 = ogr.Feature(layerDefinition)
feature2 = ogr.Feature(layerDefinition)

## Lets add the points to the feature
feature1.SetGeometry(point1)
feature2.SetGeometry(point2)
```

Notice how just now the feature and previously the layer had to get a geometry `wkbPoint` defined. After defining this geometry type, coordinates were added to the point. Then features were made based on properties of the layer. Finally we added the geometry to the feature. 

Now that we created features with geometry, we can store the features in the layer.

```{r, engine = 'python', eval=FALSE}
## Lets store the features in a layer
layer.CreateFeature(feature1)
layer.CreateFeature(feature2)

print("Extent of layer")
print(layer.GetExtent())
```

It has an extent now and it might seem we are finished now. We created geometry, features and a layer and assigned the `Esri Shapefile` driver to the layer. However the file still has to be saved. OGR doesn't have a Save() option. The shapefile is updated with the object structure when the script finished or when it is destroyed, if necessary SyncToDisk() maybe used.

```{r, engine = 'python', eval=FALSE}
## End the connection to the data source
ds.Destroy()
```

Okay nice. You created a shapefile in Python from scratch. This is just a start. Python has way more to offer for geo-science and remote sensing. Let's try some format changing and spatial processing. First try to export one of the points we created to a KML.

```{r, engine = 'python', eval=FALSE}
## Export to other formats/representations:
print("KML file export")
print(point2.ExportToKML())
```

Spatial operations are possible to. Buffer around the point and export the buffer to a GML file.

```{r, engine = 'python', eval=FALSE}
## Buffering
buffer = point2.Buffer(1,1)
print(buffer.Intersects(point1))

## More exports:
buffer.ExportToGML()
```

More spatial methods are available from the OGR module. If you are interested have a look at the methods section of this [tutorial](http://www2.geog.ucl.ac.uk/~plewis/geogg122/_build/html/Chapter4_GDAL/OGR_Python.html).

For our purposes, “ogr” will be used most of the time with vector data. Other popular drivers besides ogr are postgres (for PostgreSQL databases), KML, GML, Esri Shapefile or GeoJSON. [The python qgis cookbook](http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/loadlayer.html) gives more information on the use of the method and available drivers.

## Modify a shape file

Now let's try to add an extra point to the shape file:

```{r, engine='python', eval = FALSE}
# Drv is the ESRI shapefile driver from above and fn the filename
ds = drv.Open(fn, 1)

## Check layers and get the first layer
layernr = ds.GetLayerCount()
print(layernr)
layer = ds.GetLayerByIndex(0)
print(layer)

## Get number of features in shapefile layer
features_number = layer.GetFeatureCount()
print("number of features for this layer:", features_number)

## Get the feature definition
featureDefn = layer.GetLayerDefn()

## Create a new point
point = ogr.Geometry(ogr.wkbPoint)
point.SetPoint(0,5.6909725,50.8513682)
print(point)
## Similarly point.AddPoint(2,1)

## Create a new feature
feature = ogr.Feature(featureDefn)
feature.SetGeometry(point)

## Lets store the feature in file
layer.CreateFeature(feature)
ds.Destroy()
```

## Convert a shapefile to GeoJSON format from bash

Conversion of file formats is easier using the bash command line or even online:

See here for a web conversion example

* [https://ogre.adc4gis.com/](https://ogre.adc4gis.com/)

Using the command line
```{r, engine='bash', eval = FALSE}
ogr2ogr -f GeoJSON -t_srs crs:84 points.geojson points.shp
```

`ogr2ogr` is a command-line version of ogr and is mainly used for file conversion and/or processing. The syntax is in reverse, with new filename first and orginal filename last.

If you want to learn more about GEOJSON:

* [https://github.com/frewsxcv/python-geojson](https://github.com/frewsxcv/python-geojson)

## Visualizing spatial information

There are multiple solutions, GIS software (QGIS), web platforms (leaflet), PNG (mapnik) etc.

With leaflet a GeoJSON can be easily vizualized inside a webpage. The pythonscript below makes a html file in your working directory, which you can open in your browser: 

```{r, engine='python', eval = FALSE}
import folium
map_osm = folium.Map(location=[45.5236, -122.6750])
map_osm.save('osm.html')
```
<p>&nbsp;</p> 
<div  class="img-responsive center-block" align="center"><iframe src="osm.html" width="500" height="300" style="border: 0"></iframe></div>
<p>&nbsp;</p>

```{r, engine='python', eval = FALSE}
import folium
import os
pointsGeo=os.path.join("points.geojson")
map_points = folium.Map(location=[52,5.7],tiles='Mapbox Bright', zoom_start=6)
map_points.choropleth(geo_data=pointsGeo)
map_points.save('points.html')


```
<p>&nbsp;</p> 
<div  class="img-responsive center-block" align="center"><iframe src="points.html" width="500" height="300" style="border: 0"></iframe></div>
<p>&nbsp;</p>

Webtutorial about folium: [web-mapping-with-python-and-folium](http://pythonhow.com/web-mapping-with-python-and-folium)
 and [latest webmapping with folium](http://folium.readthedocs.io/en/latest/quickstart.html).
 
 
## Clean up
From the terminal you can remove the created "points.shp" and related files using:
```{r, engine='bash', eval=FALSE}
echo "in the terminal you can use the following bash commands to easily remove files"
echo "list all the files starting with points*"
ls points*
echo "remove all the files starting with points*"
rm -v points*
```

# What have we learned?

1. Creating or getting a writeable layer 
2. Adding fields if necessary
3. Create a feature
4. Populate the feature
5. Add the feature to the layer 
6. Close the layer
7. Modify a shapefile
8. Visualise a vector object using Python and Leaflet

# Assignment

Create your own shapefile of two point locations (e.g. a location of the GAIA building in Wageningen) in RD_New, buffer two points and make a webmap in Python. Make a well structured IPython Notebook.

- Create two points with coordinates from [Google Maps](https://www.google.nl/maps/?hl=en) or use the two points from the [R vector tutorial](https://geoscripting-wur.github.io/IntroToVector/).
- Set the projection (Lat/Long, you should know now the EPSG) and reproject to the Dutch projection RD_New
- Create a buffer around the two points of 100m
- Export the two points to shape file (that can be imported in QGIS)
- Make a webmap of your buffer with leaflet (Hint: your buffer is a polygon)

Bonus: print the distance between the two points (Hint: `Distance()`)

Upload your documented and well structured Ipython Notebook to a GitLab repository.

For making nice maps see the following [python notebooks](https://github.com/GeoScripting-WUR/PythonWeek/blob/gh-pages/Plot%20Map.ipynb).

# More information

If you have time left:

* [OGR Python tutorial](http://www2.geog.ucl.ac.uk/~plewis/geogg122/_build/html/Chapter4_GDAL/OGR_Python.html).
* [GDAL in Python](http://gdal.org/python)
* [GDAL tutorial](http://www.gdal.org/gdal_tutorial.html)
* [Geospatial Python](http://geospatialpython.com)
* [Python GDAL/OGR Cookbook](https://pcjericks.github.io/py-gdalogr-cookbook/)  

<!---

# Assignment

Upload your Python script to Github (as a .py file) uploaded within the `Python` folder part of your repository for the Python week:
Create a script that:

- defines two points in Google Earth or use the two points from the R vector excercise.
- set the projection (Lat/Long, should you know now the EPSG) and reproject to the Dutch projection (See definition in the Vector Lesson)
- draw a line between the points
- create a buffer around the points of 100m
- measure the distance between these points (i.e. the distance of that line) and print it out.
- export the two point to shape file (that can be imported in QGIS)

-->