WUR Geoscripting

Python Programming

In the previous tutorial, we learned how to set up virtual environments to write and run python code. In today’s tutorial, we will start using these environments in Positron. Next week we will work on analyzing spatial data in python.

During this tutorial, we will demonstrate different ways to visualize spatial data. This can be done using several open source packages that built upon each other. To understand the functionality of these packages and how they can be integrated, we will refer to the concepts of Object Oriented Programming (OOP). Object Oriented Programming is a way of programming where objects are fundamental building blocks that supports code modularity and reusability. OOP is a common way to re use blocks of code and add to them, especially in the development of neural networks this is a common method. Data might be loaded using a class that has some main functionality, but some custom functions can be added to make the data loader specific to a project. Therefore, before we go into the visualization part of this tutorial we will begin by explaining Object Oriented Programming.

Today’s Learning objectives

• Introduce the concept of Object Oriented Programming
• Visualize data using Matplotlib and understand the structure Matplotlib uses
• Use Cartopy to create maps
• Add geospatial data to a map using Rasterio for rasters, GeoPandas for vector

Dependencies

For the tutorial today we make use of a set of packages. Initialize a new pixi environment with the following packages in it:

pixi init python-programming
cd python-programming
pixi add python cartopy geopandas rasterio matplotlib contextily

Object-Oriented Programming in Python

Up until now in this course we have looked at R mainly as a scripting language. We call this way of programming Procedural Programming, where procedures (functions) are called in series of steps. Both Python and R can be used in another programming paradigm: Object Oriented Programming. Object Oriented Programming (OOP) is a way of programming where functionality and information is encapsulated in objects. Instead of assigning values and functions to variables individually, objects are used where both values (properties, attributes) and calculations (functions, methods) can be stored together. This offers several advantages. OOP promotes modularity and re-usability by breaking down complex problems into smaller, manageable units - the objects. These objects can be reused in various parts of the program or even in other projects, leading to more efficient, scalable and organized programming. This is especially important when working on projects containing lots of code. OOP will make your work a lot easier to understand for you and others and it will make it easier to re-use parts of the code.

How to work with objects in Python

In Python, objects are created and manipulated using classes. A class serves as a blueprint that defines the structure and behavior of an object. It brings together data (properties) and functions (methods) into a single object. To define a class in Python, we use the class keyword, followed by the name of the class. Let’s take a look at an example of a simple class called Person:

class Person:
    def __init__(self, name, age):
        self.name = name  # < this is a property
        self.age = age    # < this is also a property
    
    def greet(self):  # < This "function" is a method
        print(f"Hello, my name is {self.name} and I'm {self.age} years old.")

In Positron, you can create this class by pasting the code above into a .py file, saving it. To run this code, pick your preferred method from the previous tutorial.

However, to walk through the code interactively step by step, it’s better to use a REPL or a Jupyter Notebook. You can start a REPL or open a Jupyter Notebook, paste the code above into a cell, and run it. Note that, even though we selected the interpreter in Positron earlier for running .py files, Jupyter Notebooks run code through their own kernel, so we need to select the interpreter (python-programming in our case) again here to link the notebook’s kernel to the same environment.

Before we use the class, let’s first examine it. The Person class has two properties, name and age, as well as one method, greet. The greet method prints a greeting message that includes the person’s name and age.

In the provided code, the __init__ method is a special method known as a constructor. It is automatically called when an object (also known as an instance of a class) is created from the class. The self parameter refers to the instance of the class itself, allowing access to its properties and methods that have been assigned during instance creation or when the class was defined. Whenever a method is defined within a class (like the greet method), we give self as the first parameter.

To create an instance of the Person class, you simply call the class as if it were a function and assign the result to a variable. In REPL or Jupyter Notebook, make sure to run the cell containing the Person class first, then run the code below in later cells to create and use its instances.

person1 = Person("Alice", 25)
person2 = Person("Bob", 30)

We have created two objects, person1 and person2, which are instances of the Person class. Now we can access the properties and call the methods of these objects:

print(person1.name)  # Output: Alice
print(person2.age)   # Output: 30
person1.greet()      # Output: Hello, my name is Alice and I'm 25 years old.
person2.greet()      # Output: Hello, my name is Bob and I'm 30 years old.

Question 1: Take a look at the implementation of a geoseries object in GeoPandas. Don’t be intimidated by the amount of code! It is not necessary to understand all of it. At line 948 the plot method is defined it calls the plot_series function as defined here. What library is used for plotting and what exactly does self refer to when the plot method is defined?

Inheritance

You may have noticed that class definitions differ very slightly from what we have learned. In the GeoSeries example, the class is defined as follows:

from geopandas.base import GeoPandasBase, Series

class GeoSeries(GeoPandasBase, Series):
    pass

However, we previously learned to define a class like this:

class GeoSeries:
    pass

The difference lies in the code within the parentheses, which represents inheritance. The class will inherit all the functionality from the classes specified as arguments, in this case, GeoPandasBase and Series. The Series object refers to the pandas.Series, which contains thousands of lines of code with various functionality. Therefore, the GeoPandas GeoSeries contains all the functionality implemented in pandas, as well as those from GeoPandasBase and geopandas.Series itself. When new functionality is developed for the pandas package, this is directly available in the GeoPandas objects, since we inherit all functionality from pandas.

Phew, that’s a lot of complicated code, and it may seem overwhelming. You don’t need to understand every detail. The important thing is to grasp the value of classes, objects, and inheritance. When creating a class, we can inherit functionality from another class. How does this work in a simple example?

So, we’ve learnt that with inheritance, classes can inherit all the functionality from other classes and build upon that. Let’s create a new object called Student, which will inherit all the properties and methods from the Person class we defined earlier. In the end, a student is a person, and it possesses all its characteristics. Paste the code below to a new cell in REPL or Jupyter Notebook and run it.

class Student(Person):
    def __init__(self, name, age, student_id):
        super().__init__(name, age)
        self.student_id = student_id
        self.is_studying = False
    
    def study(self):
        self.is_studying = True
        print(f"{self.name} is studying.")

In the provided code, the Student class inherits from the Person class, which we will refer to as the superclass. By doing so, it extends the functionality of the superclass by adding a new property (student_id) and a new method (study). The super() function is used to call the superclass’s __init__ method (in this case from the Person class), allowing the subclass to initialize the inherited properties.

As a result, the Student class contains both the methods and properties inherited from the Person class, as well as the additional ones defined within the Student class:

student = Student("Eve", 22, "123456")
print(student.name)         # Output: Eve
print(student.student_id)   # Output: 123456
student.greet()             # Output: Hello, my name is Eve and I'm 22 years old.
student.study()             # Output: Eve is studying.

In this example, student is an instance of the Student class. It can access the inherited properties from the Person class, such as name, as well as the newly added property student_id and is_studying, which defaults to False. Similarly, it can invoke both the inherited method greet and the additional method study, which are specific to the Student class. The method study prints a message and sets the is_studying property to True.

Question 2: Create a new class called Teacher. This new class also inherits from Person. Define a method for the teacher that checks whether a student is studying. The student should be an input to the method.

The example of students and persons is straightforward, so let’s look at a real world example. PyTorch is a package that allows the development of neural networks. Neural networks are the networks used in deep learning, this is not the place to go in depth about deep learning, for now it is enough to understand that a neural network is a model that is trained using data, and stores information in parameters and weights. The data that is used to train the model can be text (in the case of large language models), images (in image recognition or in satellite image analysis) or other types of (spatial) data. Since we are often talking about a lot of data, a simple for loop over a path of files is often not sufficient. This is where a the definitions of a Dataset and a DataLoader come in. In this data tutorial from PyTorch, a custom dataset class is defined. This class inherits from the class Dataset, and extends its functionality with project specific data loading.

This is an often used way of working in the field of deep learning. Also using a pre-built model and extending it by adding a few layers. Your new model object would then inherit all functionality from another model object, and extend on it by adding functionality.

Visualization

Doing advanced analysis is one thing, but communicating these results is just as important. Results in raw numbers will never reach a bigger audience, and a large part of this communication depends on effective visualizations. Creating good visuals is challenging. A visual can be a graph, a colored table or an advanced infographic. In the case of geospatial data analysis, data often ends up on a map. There are many tools to visualize data, spatial or other. Maps are often made using desktop software such as QGIS, but when you want to repeatedly make similar maps, python and R can also be used. Python offers several packages to make maps, from simple static simple ones to elaborate interactive webmaps.

The most basic and one of the most used tools is Matplotlib, a general plotting package. It is used as a base for many other, more specific visualization packages. One of the core advantages of Matplotlib is that the representation of the figure is separated from the act of rendering it. This enables building increasingly sophisticated features and logic into the figure, a bit like adding many layers to a map in a Desktop GIS. Matplotlib can be used to create simple graphs but also maps. Have a look at the Python Graph Gallery to see some examples, but don’t look at the code yet! We will take you through it step by step.

Matplotlib

In the most basic form plotting is very easy. Look at the code below.

import numpy as np
from matplotlib import pyplot as plt

# Create some data
x = np.arange(-np.pi, np.pi, 0.2)
y = np.sin(x)

# Plot x against y
plt.plot(x, y)

# Show the plot
plt.show()

In this example, we are using the package NumPy to create a series of x values (values from -pi (3.1415…) to pi with steps of 0.2, check what this looks like). The y values are the sine of these values. Plotting these values is straightforward. However, as said, Matplotlib is a plotting package where a visualization object can be created before it is rendered (shown). In the example above, the rendering of the image is only done at the line plt.show(). Before this command, we can modify the plot. This allows you to add things to the visualization in steps, layering the complexity. To understand how this works, it is important to understand the hierarchy of the figure object. Have a look at the image below.

Figure 1: Matplotlib hierarchy of figure elements, source https://www.aosabook.org/en/matplotlib.html.

The basic elements are the figure and the axes objects (not to be confused with Axis objects!). The figure is like the canvas and the axes is the part of the canvas on which we will make a visualization containing for example an x-axis, y-axis (confused yet?), lines and text. Let’s build up a simple line figure as an example.

Figure 2: Matplotlib simple figure

Instead of a single plot, we can add different plots to the figure, for example two. Let’s try to add another plot with the cosine values.

# Create some data
import numpy as np 

x = np.arange(-np.pi, np.pi, 0.2)
sine = np.sin(x)
cosine = np.cos(x)

We can use the subplots method to create a figure and an array of two axes, one for each plot. Check the subplots documentation to learn more. By using the sharex and/or sharey argument, we can share axis between different plots for easy comparison between the subplots. Let’s initiate the figure and the two axes (don’t confuse axis and axes! we are initiating 2 axes, one for each plot) and let the figure share the x-axis.

import matplotlib.pyplot as plt 

# Initiate a figure with two subplots 
f, axarr = plt.subplots(2, sharex=True)

Question 3: axarr is an array. What are the elements of this array and how many elements does it consist of?

We can add a title to the figure and labels to the axes. Check this documentation to see what more you can tweak.

# Subplots are stored in an array
line = axarr[0].plot(x, sine)
axarr[0].set_title('sine plot')
axarr[1].plot(x, cosine)
axarr[1].set_title('cosine plot')
f.suptitle('This is an image of two plots', fontsize=16)

We looked now at the axes, let’s look at the axis as well. We can change the positions of the tick marks to something meaningful in the context of trigonometric functions and customize the labels with regular text or style it using LaTeX. Check if you understand what object in the hierarchy is edited and why.

# Axis label
axarr[1].set_xlabel('x')
axarr[0].set_ylabel('sin(x)')
axarr[1].set_ylabel('cos(x)')


new_ticks = np.arange(-np.pi, np.pi + 0.1, 0.25 * np.pi)
new_labels = [r"$-\pi$", r"$-\frac{3}{4}\pi$",
              r"$-\frac{1}{2}\pi$", r"$-\frac{1}{4}\pi$",
              "$0$", r"$\frac{1}{4}\pi$",
              r"$\frac{1}{2}\pi$", r"$\frac{3}{4}\pi$",
              r"$2\pi$"]
axarr[1].set_xticks(new_ticks)
axarr[1].set_xticklabels(new_labels)

Finally! our plot is done for now. We have entered a lot of commands, stacked a lot of different layers to our plot, but we cannot see it yet. using the plt.show() command we can see the result!

plt.show()

Note! We expect that you will use Positron to run this code, but for future reference, if you want to try Jupyter, some extra information. When trying to plot using a Jupyter Notebook, you may get nothing at the last step, if the cells were copy pasted step by step. This is because Jupyter automatically renders and then closes the current figure at the end of each cell, so when you run the last cell with only plt.show(), there is no active figure left to display. The easiest solution is to paste all the code into a single cell and run it together! Alternatively we can use plt.savefig('filename.png') instead of showing it. Make sure to create the plot before you run this! plt.show() closes the current plot, so calling savefig after show will result in an empty image. This is useful, otherwise we would keep adding stuff to the same plot.

two subplots with shared x-axis

One can create multiple subplots (axes) and use different plotting styles, changing for example the marker style, line style, marker size, and colors, see the example below. For the upper left subplot it is demonstrated how to add a legend; adding a label to the plotted line is essential for this.

from matplotlib import pyplot as plt

x = [1, 2, 3, 4, 5]
y = [6, 7, 8, 9, 10]

# New: define number of rows and columns of subplots and unpack them directly 
# into variables that then each contain one axes object
f, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2)

# Dashed line, label for legend, and show the legend on the subplot
ax0.plot(x, y, 'r--', label='red dashed line')
ax0.legend(loc='lower right')

# Scatter plot, using a colormap based on the y-value, changing the marker size to 35
ax1.scatter(x, y, c=y, cmap='bwr', s=35)

# Bar chart, changing the bar color to black
ax2.bar(x, y, color='k')

# Horizontal bar chart, changing the bar color to yellow
ax3.barh(x, y, color='y')
plt.show()

Matplotlib plot type examples

Question 4: In the upper right subplot, why is there no point at x=3, y=8?

There are more types of graphs available, have look at the Matplotlib documentation and play around to find out more!

Visualizing spatial data

So, plotting sine and cosines is fun and all, but this is a course about geoscripting, so how is this going to help you making maps? Well, Matplotlib really is fundamental when it comes to plotting things in python. Almost every package that can be used to create static plots (and even some dynamic ones) is based upon or uses Matplotlib, and lots of the logic you just saw will therefore be reused: the figures-, axes- and axis-objects are reused in many packages. As we saw, GeoPandas plot functionality is completely based upon Matplotlib. We will show you how to create maps using Cartopy, a geospatial wrapper around Matplotlib. Have a look at the definition of the GeoAxes. What does it inherit from? And what does this mean for its functionality?

For working with vector data we will use among others GeoPandas and for raster data we will use Rasterio packages. A more elaborate introduction to these packages will follow in the respective tutorials covering raster and vector analysis. In this tutorial we will use these packages already for reading data, if you do not entirely understand why we do certain things related to these packages, most likely this will be cleared up next week.

For this tutorial, part of the material was taken from the Project Pythia website, an excellent source for more information and tutorials about working with python!

Cartopy

As said, Cartopy is basically adding the spatial component to Matplotlib. Therefore, a lot of the logic will be familiar. By adding spatial information (a coordinate reference system) to an Axes (turning it into a GeoAxes) we can reference our spatial data to each other. Additionally, Cartopy has a built-in module to handle the referencing cartopy.crs. Let’s import these libraries and modules.

import matplotlib.pyplot as plt
from cartopy import crs as ccrs
from cartopy import feature as cfeature

Let’s start by creating a map of the world. We do this by generating 1 subplot (so just a plot) with the Plate Carrée projection, a projection where every point is spaced out equally in terms of degrees.

fig = plt.figure(figsize=(11, 8.5))
ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=-75))
ax.set_title("A Geo-referenced subplot, Plate Carree projection")

We don’t see anything yet! That’s because we have not put anything on the map. Let’s add the coastline for some spatial context, which can be done by calling a method of the GeoAxes object ax.coastlines.

ax.coastlines()

coastlines is a special case, apparently showing the coastlines happened so often that a special function was defined. Not everything is this simple sadly… In the cartopy documentation we can see that there exist pre defined features that we can add using the add_feature method from a GeoAxes. ax.add_feature(cfeature.COASTLINE, linewidth=0.3, edgecolor='black') does the same as ax.coastlines() effectively (don’t believe it? Check the source! ). Have a look and play around with adding other features! Using the linewidth andedgecolor arguments we can style the map. Don’t forget to plt.show() to see the map!

Question 5: Create a worldwide map with 3 different features, each styled differently. Also add the stockimage to the map. Use the documentation to see how.

We have now step by step built up a map in a Plate Carrée projection system. However, this projection has some major issues, have you seen how big Antarctica would be at the equator?! Let’s quickly create another map in another projection. In the documentation we can find a large list of projections that can be used.

fig = plt.figure(figsize=(11, 8.5))
projLae = ccrs.LambertAzimuthalEqualArea(central_longitude=0.0, central_latitude=0.0)
ax = plt.subplot(1, 1, 1, projection=projLae)
ax.set_title("Lambert Azimuthal Equal Area Projection")
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue')
plt.show()

Local maps

We have seen how to make worldwide maps, let’s now make a map, closer to home, with our own data. The polygon data that is shown here is read by GeoPandas, directly from a url. The add_geometries method reads the geometries from a GeoDataFrame, but can also read geometries from Shapely (more about this in the Python Vector tutorial). The set_extent method is used to define an area of interest, basically it sets the top and bottom left and right corners, so that only the area of interest is shown. Have a look at the code below, the comments explain more line by line.

import geopandas as gpd

gdf = gpd.read_file('https://raw.githubusercontent.com/GeoScripting-WUR/PythonProgramming/master/data/gadm41_NLD_2.json')

# Using the dutch coordinate reference system RDNew (epsg code 28992)
crs_rd = ccrs.epsg(28992)

# The data is in another projection than our plot, reprojection to RDNew
gdf = gdf.to_crs(28992)

# Initiate the plot, a little bigger than before
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, projection=crs_rd)
ax.set_title('The municipalities of NL')

# Draw gridlines
gl = ax.gridlines(
    draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--'
)

# Set the extent to the extent of the municipalities
min_x, max_x, min_y, max_y = gdf.total_bounds
ax.set_extent((min_x, max_x, min_y, max_y), crs=crs_rd)

# ax.set_extent(gdf.total_bounds, crs=crs_rd) would do this in one step,
# but the coordinates can be defined separately as well in this order!

# Add the geometries to the map
ax.add_geometries(gdf["geometry"], crs=crs_rd, edgecolor = 'black', facecolor = 'None')

plt.show()

Basemaps

Our map of the municipalities is looking a bit empty though, an empty white shape does not give much context to work with. It would be nice to see what’s actually there: roads, cities, water, the kind of context you get “for free” on Google Maps. This is called a basemap.

Basemaps are often supplied using map tiles. Cartopy can fetch map tiles itself through cartopy.io.img_tiles, but it is fiddly to work with, you need to pick zoom levels by hand and several of the free tile providers it relies on have been discontinued over the years. A more lightweight, purpose built package for this is contextily. Remember that a GeoAxes inherits from a regular Matplotlib Axes? That means any tool that works on a normal Axes, like contextily, works on our GeoAxes too. So, we can take the exact map we just built above and simply add a basemap to it.

import contextily as cx

# add_basemap fetches tiles for the current extent of ax, and reprojects
# them on the fly to whatever crs we give it
cx.add_basemap(ax, crs=28992, zorder=-1)

plt.show()

Keep in mind that in this code, we add to the already existing ax, reusing the earlier plot from the municipalities. This above codeblock relies on the fact that you ran that code before.

The image tiles that are used as a basemap in this example are by default available in the web mercator projection (epsg:3857), the crs argument allows for on the fly reprojection to whatever we need, in this case the dutch RD New coordinate system.

The zorder argument controls the stacking order of the different layers on our GeoAxes, just like layers in a GIS. Giving the basemap a low zorder makes sure it stays behind the municipality polygons we added earlier, instead of covering them up.

Question 6: By default, contextily uses OpenStreetMap tiles. Have a look at the contextily documentation to find another tile provider, and use it instead.

Plotting rasters

In the case of rasters, we will make use in another way of the widespread use of Matplotlib and the GeoAxes objects. For rasters we will make use of the package Rasterio to handle raster files. In this package a plot module is defined, that allows for the plotting of rasters. Instead of defining an axes and adding the raster as a feature to the plot, we will create the Figure and GeoAxes objects using Cartopy, and pass them to Rasterio’s plot.show function. Other features that we might want to add, we can add still to the same GeoAxes, in that way everything will be added to the same canvas. Have a look at the code below.

import requests
import io 
import zipfile 
import rasterio
from rasterio.plot import show
import numpy as np


# These first 4 lines download and unzip a landsat8 image 
# It is not necessary to understand these lines. 
url = 'https://github.com/GeoScripting-WUR/VectorRaster/releases/download/tutorial-data/landsat8.zip'
resp = requests.get(url)
zf = zipfile.ZipFile(io.BytesIO(resp.content))
zf.extractall('./')

#The landsat image is projected in UTM31N, let's use that projection
crs_utm = ccrs.epsg(32631)

# Initiate the area of interest
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, projection=crs_utm)
ax.set_title('The municipalities of NL')

# Read the GeoJSON with GeoPandas.
gdf = gpd.read_file('https://raw.githubusercontent.com/GeoScripting-WUR/PythonProgramming/master/data/gadm41_NLD_2.json')
# This time reproject it to UTM31
gdf = gdf.to_crs(32631)
gdf.plot(ax=ax,edgecolor='white', color = 'None')

# Open the raster using Rasterio, more about Rasterio next week!
dataset = rasterio.open('./LC81970242014109LGN00.tif')

# This is a Landsat 8 image with 7 bands. For a true-color image we want
# the Red, Green and Blue bands, which for Landsat 8 are bands 4, 3 and 2.
rgb = dataset.read((4, 3, 2))

# Most valid pixels fall between 0 and 1500, with a long tail. Clipping to 
# that range and rescaling to 0-1 to make the image brighter
rgb = np.clip(rgb, 0, 1500) / 1500


show(rgb, transform=dataset.transform, ax=ax)

Remember to show or save the image! For this, run plt.show() to show the image, or plt.savefig('filename.png') to save the image.

Adding legends

We have seen now how to add spatial data to a plot and how to add context with a basemap. Finally, putting data on a map is one thing, but to make the maps interpretable we need to tell the reader what a color means. For this we will need to add a legend to the map. Cartopy itself does not allow us to do that, but the friendly folks from GeoPandas have made this easy for us when plotting vector data, and we will reuse some Matplotlib functionality when plotting raster data.

Legends on vector data.

Before showing how to add a legend, have a look at the GeoDataFrame.plot documentation. There are quite some different arguments that we can pass to the plot method of a GeoDataFrame object. Most importantly, there is the legend argument that lets GeoPandas add a legend to a plot. See the codeblock below to see how we add a legend after plotting a GeoDataFrame with polygons.

# Vector plotting with legend
fig2 = plt.figure(figsize=(15, 15))
ax2 = plt.subplot(1, 1, 1, projection=crs_utm)
ax2.set_title('Municipalities of NL, colored by province')

# GeoDataFrame.plot() works directly on a GeoAxes when the data's CRS
# matches the axes' projection, and (unlike add_geometries) it builds a
# legend for us when we give it a categorical column
gdf.plot(
    ax=ax2,
    column='NAME_1',
    categorical=True,
    legend=True, # This is where we add the legend
    edgecolor='black',
    legend_kwds={'loc': 'lower right', 'title': 'Province', 'fontsize': 8}, # Here we add some configuration to the for the placing and fontsize in the legend 
)
ax2.set_extent((min_x, max_x, min_y, max_y), crs=crs_utm)

plt.show()

There are some caveats and for a long time there were issues and a wishlist with features to change this behaviour. For example, in the current version it is cumbersome to add legend entries for other geometry types than polygons. The legend will always show the color rather than the geometry type (including line thickness or point styling). In the next GeoPandas version however, this will be fixed, see this merged pull request. It contains quite some new functionality, it is also a good illustration of how open source software development works!

Legends on raster data

For adding a legend, for example a colorbar, to a raster plot there is no built-in functionality in Rasterio. Instead, we will use Matplotlib’s colorbar functionality. For the example, see below, we cannot use the full color RGB image from above, but we need to display a single band. We will show the Digital Number (DN) from the nir band in the legend. For the plotting we will use Rasterio’s show, as explained above.

from rasterio.plot import show


# Acquiring the data (ignore this code)
url = 'https://github.com/GeoScripting-WUR/VectorRaster/releases/download/tutorial-data/landsat8.zip'
resp = requests.get(url)
zf = zipfile.ZipFile(io.BytesIO(resp.content))
zf.extractall('./')
dataset = rasterio.open('./LC81970242014109LGN00.tif')

# using the same CRS as the data 
crs_utm = ccrs.epsg(32631)

# Initializing the figure
fig3 = plt.figure(figsize=(15, 15))

# Creating a geoaxes and setting its title
ax3 = plt.subplot(1, 1, 1, projection=crs_utm)
ax3.set_title('Landsat 8 - near-infrared band (band 5)')


# Add a single band to the plot, so each value maps to one color - 
# this is what a colorbar visualizes. 
# An RGB composite has no single value to map.
nir = dataset.read(5)

# use Rasterio's show to add the nir band to ax3
show(nir, transform=dataset.transform, ax=ax3, cmap='viridis')

# show() returns the ax, not the image object needed for a colorbar.
# That image is stored on the ax though, so we grab it from there.
im = ax3.images[0]
fig3.colorbar(im, ax=ax3, label='Digital Number (DN)', shrink=0.7)

plt.show()

What have we learned?

You finished this tutorial, well done! We started with a introduction about object oriented programming and you now know what objects are, how they are implemented in python and how they can inherit functionality from each other. In this way we can stand on top of the shoulders of giants, we do not have to write the same code somebody else already has.

Next Matplotlib was introduced, you now know what the structure of a Matplotlib plot is, how different elements are organized on a figure and how to add data to a plot.

Building on that, we looked at Cartopy, building on top of Matplotlib, illustrating object oriented programming and showing its power. We made maps using all built in functionality, and we also saw how to add our own data to maps both vector and raster. Because a GeoAxes is still a Matplotlib Axes, we could also drop in contextily to add a basemap without any extra work, another example of reusing code somebody else already wrote. We used Cartopy’s GeoAxes and passed them to both GeoPandas and Rasterio’s plotting functionality to add layers to the plots. We also saw how to add legends to the plots. Keep an eye out to the next GeoPandas version for improved legend functionality!

What you learned today is only a tip of the iceberg. For more elaborate plots and other examples, visit the Pythia project and the Cartopy gallery, both listed below!

More info

• Official Python tutorial
• Python Style guide
• Python 3 Cheatsheet
• Overview Python package Cheatsheets
• Project Pythia
• Cartopy examples