Week 3: Python for geo-scripting

WUR Geoscripting WUR logo

Vector data handling with Python

Jan Verbesselt, Jorge Mendes de Jesus, Aldo Bergsma, Dainius Masiliunas - 22 January, 2018

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. 
python
try:
  from osgeo import ogr
except:
  import ogr
From the terminal you can also type the following to obtain an overview of supported formats: 
echo "Terminal command to obtain an overview of supported formats"
ogrinfo --formats

More info about the different OGR formats.

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.

## Display all your conda environments
conda info --envs

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

## 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.

## 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.

from osgeo import ogr
# Create a point geometry
wkt = "POINT (173914.00 441864.00)"
pt = ogr.CreateGeometryFromWkt(wkt)
print(pt)
# print(help(osgeo.ogr))
## POINT (173914 441864)

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

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.

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

from osgeo import osr
## Spatial reference
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(4326)  # from EPSG - Lat/long

Reproject a point

## 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.

This can be visualised as a pyramid of elements:

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
## 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.

## 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.

# 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.

## 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)

## 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.

## 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.

## 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.

## 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.

## 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.

## 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.

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 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:

# 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

Using the command line 
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:

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:

import folium
map_osm = folium.Map(location=[45.5236, -122.6750])
map_osm.save('osm.html')

 

  

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')

 

 

Webtutorial about folium: web-mapping-with-python-and-folium and latest webmapping with folium.

Clean up

From the terminal you can remove the created "points.shp" and related files using: 
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 or use the two points from the R vector tutorial.
  • 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.

More information

If you have time left: