diff --git a/.gitignore b/.gitignore index d43d1c4..d371004 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ # pixi environments .pixi/* !.pixi/config.toml +__pycache__ \ No newline at end of file diff --git a/PythonProgrammig.Rproj b/PythonProgrammig.Rproj deleted file mode 100644 index 8e3c2eb..0000000 --- a/PythonProgrammig.Rproj +++ /dev/null @@ -1,13 +0,0 @@ -Version: 1.0 - -RestoreWorkspace: Default -SaveWorkspace: Default -AlwaysSaveHistory: Default - -EnableCodeIndexing: Yes -UseSpacesForTab: Yes -NumSpacesForTab: 2 -Encoding: UTF-8 - -RnwWeave: Sweave -LaTeX: pdfLaTeX diff --git a/TestGround/cartopy.py b/TestGround/cartopy.py index 376e52b..5aca7b1 100644 --- a/TestGround/cartopy.py +++ b/TestGround/cartopy.py @@ -1,39 +1,115 @@ +import io +import zipfile + import geopandas as gpd +import matplotlib.pyplot as plt +import pandas as pd +import requests +import rasterio +from cartopy import crs as ccrs +from rasterio.plot import show +import contextily as cx -gdf = gpd.read_file('https://raw.githubusercontent.com/GeoScripting-WUR/PythonProgramming/master/data/gadm41_NLD_2.json') +# gdf = gpd.read_file('https://raw.githubusercontent.com/GeoScripting-WUR/PythonProgramming/master/data/gadm41_NLD_2.json') -# Using the dutch coodinate reference system RDNew (epsg code 28992) -crs = ccrs.epsg(28992) +# # Using the dutch coodinate reference system RDNew (epsg code 28992) +# crs = ccrs.epsg(28992) -# The data is in another projection as our plot, reprojection to RDnew -gdf = gdf.to_crs(28992) +# # The data is in another projection as our plot, reprojection to RDnew +# gdf = gdf.to_crs(28992) -# Initiate the plot, a little bigger then before -fig = plt.figure(figsize=(15, 15)) -ax = plt.subplot(1, 1, 1, projection=crs) -ax.set_title('The municipalities of NL') +# # Add an attribute with the area of each municipality. RDNew (28992) is a +# # projected CRS in meters, so .area gives square meters straight away +# gdf['area_km2'] = gdf.geometry.area / 1_000_000 -# Draw gridlines -gl = ax.gridlines( - draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--' -) +# # Categorize the areas into 10 equal-width classes, one per municipality +# gdf['area_class'] = pd.cut(gdf['area_km2'], bins=10) + +# # Initiate the plot, a little bigger then before +# fig = plt.figure(figsize=(15, 15)) +# ax = plt.subplot(1, 1, 1, projection=crs) +# 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) -# 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) +# # ax.set_extent(gdf.total_bounds, crs=crs) would do this in one step, +# # but the coordinates can be defined seperately as well in this order! -# ax.set_extent(gdf.total_bounds, crs=crs) would do this in one step, -# but the coordinates can be defined seperately as well in this order! +# # Color each municipality by its area class. add_geometries can't do this +# # (one edgecolor/facecolor for the whole collection), so we use +# # GeoDataFrame.plot() instead - it works directly on a GeoAxes when the +# # data's CRS matches the axes' projection, and it builds a legend for us +# # when given a categorical column +# gdf.plot( +# ax=ax, +# column='area_class', +# categorical=True, +# legend=True, +# cmap='viridis', +# edgecolor='black', +# legend_kwds={'loc': 'lower right', 'title': 'Municipality area (km²)', 'fontsize': 8}, +# ) -# Add the geometries to the map -ax.add_geometries(gdf["geometry"], crs=crs, edgecolor = 'black', facecolor = 'None') +# # 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=crs, zorder=-1) + +# plt.show() + + +# ---- Vector legend example (categorical) ---- + +fig2 = plt.figure(figsize=(15, 15)) +ax2 = plt.subplot(1, 1, 1, projection=crs) +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, + edgecolor='black', + legend_kwds={'loc': 'lower right', 'title': 'Province', 'fontsize': 8}, +) +ax2.set_extent((min_x, max_x, min_y, max_y), crs=crs) plt.show() -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=crs, zorder=-1) +# ---- Raster legend example (continuous colorbar) ---- + +url = 'https://github.com/GeoScripting-WUR/VectorRaster/releases/download/tutorial-data/landsat8.zip' +resp = requests.get(url, "data.zip") +zf = zipfile.ZipFile(io.BytesIO(resp.content)) +zf.extractall('./') + +crsUtm = ccrs.epsg(32631) + +fig3 = plt.figure(figsize=(15, 15)) +ax3 = plt.subplot(1, 1, 1, projection=crsUtm) +ax3.set_title('Landsat 8 - near-infrared band (band 5)') -plt.show() \ No newline at end of file +dataset = rasterio.open('./LC81970242014109LGN00.tif') + +# A single band, 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) + +show(nir, transform=dataset.transform, ax=ax3, cmap='viridis') + +# show() returns the ax, not the image mappable, but the AxesImage it +# just drew is still available on the ax afterwards +im = ax3.images[0] +fig3.colorbar(im, ax=ax3, label='Digital Number (DN)', shrink=0.7) + +plt.show() diff --git a/_quarto.yml b/_quarto.yml new file mode 100644 index 0000000..de29693 --- /dev/null +++ b/_quarto.yml @@ -0,0 +1,7 @@ +# _quarto.yml +project: + type: website + output-dir: _site + +website: + title: "R & Python Basics" \ No newline at end of file diff --git a/index.Rmd b/index.Rmd deleted file mode 100644 index 41f15bc..0000000 --- a/index.Rmd +++ /dev/null @@ -1,452 +0,0 @@ ---- -pagetitle: "Tutorial 9: Python Programing" -author: "Arno Timmer, Jan Verbesselt, Jorge Mendes de Jesus, Aldo Bergsma, Johannes Eberenz, Dainius Masiliunas, David Swinkels, Judith Verstegen, Corné Vreugdenhil" -date: "`r format(Sys.time(), '%d %B, %Y')`" -output: - rmdformats::html_clean: - title: "Tutorial 9: Python Programing" - theme: "simplex" - highlight: zenburn - menu: FALSE - theme.chooser: TRUE - highlight.chooser: TRUE ---- - -```{css, echo=FALSE} -@import url("https://netdna.bootstrapcdn.com/bootswatch/3.0.0/simplex/bootstrap.min.css"); -.main-container {max-width: none;} -div.figcaption {display: none;} -pre {color: inherit; background-color: inherit;} -code[class^="sourceCode"]::before { - content: attr(class); - display: block; - text-align: right; - font-size: 70%; -} -code[class^="sourceCode r"]::before { content: "R Source";} -code[class^="sourceCode python"]::before { content: "Python Source"; } -code[class^="sourceCode bash"]::before { content: "Bash Source"; } -``` - -[WUR Geoscripting](https://geoscripting-wur.github.io/) WUR logo - -# Python Programing - -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 VS Code, which will serve as our default IDE. Next week we will work with spatial data analysis which more often than not results in data being displayed on a map. - -During this tutorial, we will demonstrate different ways to visualize data on a map in both static and interactive formats. 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. Since OOP is commonly used in the development of open source packages, we can reuse and adapt already developed code and its functionality. Therefore, before we go into the visualization part of this tutorial we will begin by explaining Object Oriented Programming. - - -## Today’s Learning objectives - -- Vizualize data using Matplotlib and understand the structure matplotlib uses -- Use cartopy's built in functionality 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. Create a new yaml file (`env.yaml` for example) and paste the following environment definition in it. -``` -name: python-programming -dependencies: - - cartopy - - python - - spyder - - geopandas - - rasterio - - matplotlib -``` - -In the terminal of VS Code, create the environment using `mamba env create --file env.yaml`. Then activate this environment using `source activate python-programming`. Then select `python-programming` as interpreter, just as we did in the previous tutorial. - - -# 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`: - -```{python, eval=FALSE} -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 VS Code, you can create this class by pasting the code above into a .py file, saving it, and then running it in the terminal with `python3 your_class.py`. - -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 VS Code 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 instance of an 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 `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. - -```{python, eval=FALSE} -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: - -```{python, eval=FALSE} -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. -``` - -This example is straightforward, but keep in mind that classes can become more complex. - -```{block, type="alert alert-success"} -> **Question 2**: Take a look at the [implementation of a geoseries](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/geoseries.py#L79) 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](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/geoseries.py#L948) it calls the [`plot_series` function as defined here](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/plotting.py#L313). What library is used for plotting and what exactly does `self` refer to when the `plot` method is defined? -``` - -## Inheritence - -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: - -```{python, eval=FALSE} -class GeoSeries(GeoPandasBase, Series): -``` - -However, we previously learned to define a class like this: - -```{python, eval=FALSE} -class GeoSeries: -``` - -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](https://github.com/pandas-dev/pandas/blob/main/pandas/core/series.py#L243C15-L243C15), 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. - -```{python, eval=FALSE} -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 `Pesron` 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: - -```{python, eval=FALSE} -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`. - -```{block, type="alert alert-success"} -> **Question 3**: 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. -``` - -# Visualization - -Communicating research results is challenging without good visualizations like graphs or more elaborate infographics, and in the case of geospatial data analysis, a map. There are many tools to vizualize data using python, some of them can become very elaborate. -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 tailored 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 GIS. Matplotlib can be used to create simple graphs but also maps. Have a look at the [Python Graph Gallery](https://python-graph-gallery.com/matplotlib/) 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 (again, you can try running this in both REPL or Jupyter Notebook, or your own preferred IDEs): - -```{Python,engine.path='/usr/bin/python3'} -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.14) 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 vizualization 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 to add things to the vizualization 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. - -
- - - -
Matplotlib hierarchy of figure elements, source .
- -
- -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, lines and text. Let's build up a simple line figure as an example. - -matplotlib simple figure - -Note that the behavior of `plt.show()` depends on how you run the script. If you are using Spyder and want `plt.show()` to work: - -1. Go to Tools. -2. Go to Preferences. -3. Select IPython console. -4. Go to Graphics tab. -5. In the Graphics backend section, select Automatic as the backend type. -6. Restart your kernel (go to Console > Restart kernel). - -Instead of a single plot, we can add different plots to the figure, for example two. Let's try to add another pot with the cosine values. - -```{Python,engine.path='/usr/bin/python3'} -# Create some data -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](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html) 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. - -```{Python,engine.path='/usr/bin/python3'} -# Initiate a figure with two subplots -f, axarr = plt.subplots(2, sharex=True) - -``` - -```{block, type="alert alert-success"} -> **Question 4**: `axarr` is an array. What are the elements of this array and how many elements does is exist out of? -``` - -We can add a title to the figure and labels to the axes. Check [this documentation](https://matplotlib.org/stable/api/axes_api.html) to see what more you can tweak. - -```{Python,engine.path='/usr/bin/python3'} -# 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 a 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. - -```{Python,engine.path='/usr/bin/python3'} -# 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! - -```{Python,engine.path='/usr/bin/python3'} -plt.show() -``` - - -Note that if you paste the above code blocks cell by cell in a REPL, you will see the final figure as expected. In Jupyter Notebook, however, you may get nothing at the last 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 and 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. - -```{Python,engine.path='/usr/bin/python3'} -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 - -```{block, type="alert alert-success"} -> **Question 5**: 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](https://matplotlib.org/stable/plot_types/index.html) and play around to find out more! - - -## Vizualizing spatial data - -So, plotting sine and cosines is fun and all, but this is a course about **geo**scripting, 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`](https://github.com/SciTools/cartopy/blob/main/lib/cartopy/mpl/geoaxes.py#L354). 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](https://foundations.projectpythia.org/core/cartopy/cartopy.html) 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. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -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](https://en.wikipedia.org/wiki/Equirectangular_projection) projection, a projection where every point is spaced out equally in terms of degrees. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -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`. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -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](https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html) 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()` (don't believe it? [Check the source](https://github.com/SciTools/cartopy/blob/75939f9e81ac67838c52ce80230e4431badbeace/lib/cartopy/mpl/geoaxes.py#L609)! ) effectively. Have a look and play around with adding other features! Using the `linewidth` and`edgecolor` arguments we can style the map. Don't forget to `plt.show()` to see the map! - -```{block, type="alert alert-success"} -> **Question 6**: Create a worldwide map with 3 different features, each styled differently. Also add the stockimage to the map. Use [the documentation](https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/geoaxes.html?highlight=stock#cartopy.mpl.geoaxes.GeoAxes.stock_img) to see how. -``` - -We have now step by step built up a map in a Pate Carree projection system. However, this projection has some major issues, have you seen [how big Antartica is at the equator](https://www.thetruesize.com)?! Let's quickly create another map in another projection. In the [documentation](https://scitools.org.uk/cartopy/docs/latest/reference/projections.html) we can find a large list of projections that can be used. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -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() -``` - -## Smaller maps -We have seen how to make worldwide maps, let's now make a smaller map 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 Shapelt (more about this in 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. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -import geopandas as gpd - -gdf = gpd.read_file('https://raw.githubusercontent.com/GeoScripting-WUR/PythonProgramming/master/data/gadm41_NLD_2.json') - -# Using the dutch coodinate reference system RDNew (epsg code 28992) -crs = ccrs.epsg(28992) - -# The data is in another projection as our plot, reprojection to RDnew -gdf = gdf.to_crs(28992) - -# Initiate the plot, a little bigger then before -fig = plt.figure(figsize=(15, 15)) -ax = plt.subplot(1, 1, 1, projection=crs) -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) - -# ax.set_extent(gdf.total_bounds, crs=crs) would do this in one step, -# but the coordinates can be defined seperately as well in this order! - -# Add the geometries to the map -ax.add_geometries(gdf["geometry"], crs=crs, edgecolor = 'black', facecolor = '#FFFFFF') - -plt.show() -``` - -## 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](https://github.com/rasterio/rasterio/blob/7c83de410b7b812e6552e4b76ce236bb6213c80d/rasterio/plot.py#L34) 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` `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. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -import requests -import io -import zipfile -import rasterio -from rasterio.plot import show - -# 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, "data.zip") -zf = zipfile.ZipFile(io.BytesIO(resp.content)) -zf.extractall('./') - -#The landsat image is projected in UTM31N, let's use that projection -crs = ccrs.epsg(32631) - -# Initiate the area of interest -fig = plt.figure(figsize=(15, 15)) -ax = plt.subplot(1, 1, 1, projection=crs) -ax.set_title('The municipalities of NL') - -# Read the GeoJson again 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') -show(dataset, ax=ax, cmap='gist_ncar') -``` - -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. - -# What have we learned? - -You finished this tutorial well done! -We started with a short 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 it's power. We made basic maps using all built in functionality, and we also saw how to add our own data to maps both vector and raster. - -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](https://docs.Python.org/3/contents.html) -- [Python Style guide](https://www.python.org/dev/peps/pep-0008/) -- [Python 3 Cheatsheet](https://ugoproto.github.io/ugo_py_doc/py_cs/) -- [Overview Python package Cheatsheets](https://www.datacamp.com/community/data-science-cheatsheets?tag=python) -- [Project Pythia](https://foundations.projectpythia.org/core/cartopy/cartopy.html) -- [Cartopy examples](https://scitools.org.uk/cartopy/docs/latest/gallery/index.html) diff --git a/index.html b/index.html deleted file mode 100644 index 9cb458d..0000000 --- a/index.html +++ /dev/null @@ -1,979 +0,0 @@ - - - - - - - - - - - -Tutorial 9: Python Programing - - - - - - - - - - - - - - - - - - - - - -
- - -
- -
-
-
- - - -
- -
-
Author
-
-

Arno Timmer, Jan Verbesselt, Jorge Mendes de Jesus, Aldo Bergsma, Johannes Eberenz, Dainius Masiliunas, David Swinkels, Judith Verstegen, Corné Vreugdenhil

-
-
- -
-
Published
-
-

August 21, 2026

-
-
- - -
- - - -
- - -

WUR Geoscripting

-
-

Python Programing

-

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 in both static and interactive formats. 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, especialy in the development of neural netowrks this a common method. Data might be loadeed 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

-
    -
  • Vizualize data using Matplotlib and understand the structure matplotlib uses
  • -
  • Use cartopy’s built in functionality 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 VS Code 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 instance of an 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 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 2: 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?

-
-
-
-
-

Inheritence

-

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.Seriesitself. 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 Pesron 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 3: 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 sattelite 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 it’s 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

-

Communicating research results is challenging without good visualizations like graphs or more elaborate infographics, and in the case of geospatial data analysis, a map. There are many tools to vizualize data using python, some of them can become very elaborate.

-

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 tailored 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 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.14) 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 vizualization 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 to add things to the vizualization 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, lines and text. Let’s build up a simple line figure as an example.

-
-
-

-
matplotlib simple figure
-
-
-

Note that the behavior of plt.show() depends on how you run the script. If you are using Spyder and want plt.show() to work:

-
    -
  1. Go to Tools.
  2. -
  3. Go to Preferences.
  4. -
  5. Select IPython console.
  6. -
  7. Go to Graphics tab.
  8. -
  9. In the Graphics backend section, select Automatic as the backend type.
  10. -
  11. Restart your kernel (go to Console > Restart kernel).
  12. -
-

Instead of a single plot, we can add different plots to the figure, for example two. Let’s try to add another pot 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 4: axarr is an array. What are the elements of this array and how many elements does is exist out 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 a 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 that if you paste the above code blocks cell by cell in a REPL, you will see the final figure as expected. In Jupyter Notebook, however, you may get nothing at the last 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 and 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 5: 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!

-
-
-

Vizualizing 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() (don’t believe it? Check the source! ) effectively. 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 6: 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 Pate Carree projection system. However, this projection has some major issues, have you seen how big Antartica is 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()
-
-
-
-

Smaller maps

-

We have seen how to make worldwide maps, let’s now make a smaller map 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 Shapelt (more about this in 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 coodinate reference system RDNew (epsg code 28992)
-crs = ccrs.epsg(28992)
-
-# The data is in another projection as our plot, reprojection to RDnew
-gdf = gdf.to_crs(28992)
-
-# Initiate the plot, a little bigger then before
-fig = plt.figure(figsize=(15, 15))
-ax = plt.subplot(1, 1, 1, projection=crs)
-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)
-
-# ax.set_extent(gdf.total_bounds, crs=crs) would do this in one step, 
-# but the coordinates can be defined seperately as well in this order! 
-
-# Add the geometries to the map
-ax.add_geometries(gdf["geometry"], crs=crs, edgecolor = 'black', facecolor = '#FFFFFF')
-
-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.

-

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=crs, zorder=-1)
-
-plt.show()
-
-

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 7: 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 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, "data.zip")
-zf = zipfile.ZipFile(io.BytesIO(resp.content))
-zf.extractall('./')
-
-#The landsat image is projected in UTM31N, let's use that projection
-crs = ccrs.epsg(32631)
-
-# Initiate the area of interest
-fig = plt.figure(figsize=(15, 15))
-ax = plt.subplot(1, 1, 1, projection=crs)
-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.

-
-
-
-

What have we learned?

-

You finished this tutorial, well done! We started with a short 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 it’s power. We made basic 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.

-

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

- -
- -
- - -
- - - - - \ No newline at end of file diff --git a/index.qmd b/index.qmd index 2dd19ea..c837e55 100644 --- a/index.qmd +++ b/index.qmd @@ -1,5 +1,5 @@ --- -pagetitle: "Tutorial 9: Python Programing" +pagetitle: "Tutorial 9: Python Programming" author: "Arno Timmer, Jan Verbesselt, Jorge Mendes de Jesus, Aldo Bergsma, Johannes Eberenz, Dainius Masiliunas, David Swinkels, Judith Verstegen, Corné Vreugdenhil" date: today format: @@ -9,6 +9,7 @@ format: toc: true toc-location: left css: styles.css + lightbox: true execute: eval: false echo: true @@ -16,18 +17,19 @@ execute: [[WUR Geoscripting](https://geoscripting-wur.github.io/)]{.page-header-title} -# Python Programing +# 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 in both static and interactive formats. 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, especialy in the development of neural netowrks this a common method. Data might be loadeed 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. +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 -- Vizualize data using Matplotlib and understand the structure matplotlib uses -- Use cartopy's built in functionality to create maps -- Add geospatial data to a map using rasterio for rasters, geopandas for vector +- 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: @@ -54,13 +56,13 @@ class Person: 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. +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 VS Code 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. +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 instance of an 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 `greet` method), we give `self` as the first parameter. +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. @@ -79,10 +81,10 @@ person2.greet() # Output: Hello, my name is Bob and I'm 30 years old. ``` ::: {.alert .alert-success} -> **Question 2**: Take a look at the [implementation of a geoseries](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/geoseries.py#L79) 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](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/geoseries.py#L948) it calls the [`plot_series` function as defined here](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/plotting.py#L313). What library is used for plotting and what exactly does `self` refer to when the `plot` method is defined? +> **Question 1**: Take a look at the [implementation of a geoseries](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/geoseries.py#L79) 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](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/geoseries.py#L948) it calls the [`plot_series` function as defined here](https://github.com/geopandas/geopandas/blob/80edc868454d3fae943b734ed1719c2197806815/geopandas/plotting.py#L313). What library is used for plotting and what exactly does `self` refer to when the `plot` method is defined? ::: -## Inheritence +## 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: @@ -100,7 +102,7 @@ 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](https://github.com/pandas-dev/pandas/blob/main/pandas/core/series.py#L243C15-L243C15), 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. +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](https://github.com/pandas-dev/pandas/blob/main/pandas/core/series.py#L243C15-L243C15), 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? @@ -118,7 +120,7 @@ class Student(Person): 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 `Pesron` class), allowing the subclass to initialize the inherited properties. +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: @@ -133,19 +135,19 @@ 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`. ::: {.alert .alert-success} -> **Question 3**: 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. +> **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](https://docs.pytorch.org) 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 sattelite 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](https://docs.pytorch.org/tutorials/beginner/basics/data_tutorial.html#creating-a-custom-dataset-for-your-files) from PyTorch, a custom dataset class is defined. This class inherits from the class [`Dataset`](https://docs.pytorch.org/docs/2.13/data.html#torch.utils.data.Dataset), and extends it's functionality with project specific data loading. +The example of students and persons is straightforward, so let's look at a real world example. [PyTorch](https://docs.pytorch.org) 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](https://docs.pytorch.org/tutorials/beginner/basics/data_tutorial.html#creating-a-custom-dataset-for-your-files) from PyTorch, a custom dataset class is defined. This class inherits from the class [`Dataset`](https://docs.pytorch.org/docs/2.13/data.html#torch.utils.data.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 -Communicating research results is challenging without good visualizations like graphs or more elaborate infographics, and in the case of geospatial data analysis, a map. There are many tools to vizualize data using python, some of them can become very elaborate. +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 tailored 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 GIS. Matplotlib can be used to create simple graphs but also maps. Have a look at the [Python Graph Gallery](https://python-graph-gallery.com/matplotlib/) to see some examples, but don't look at the code yet! We will take you through it step by step. +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](https://python-graph-gallery.com/matplotlib/) to see some examples, but don't look at the code yet! We will take you through it step by step. ## Matplotlib @@ -166,7 +168,7 @@ plt.plot(x, y) plt.show() ``` -In this example, we are using the package numpy to create a series of x values (values from -pi (-3.14) 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 vizualization 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 to add things to the vizualization 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. +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. ::: {#fig-mpl-structure layout-ncol="2"} ![](images/mpl_structure.png) @@ -176,21 +178,11 @@ In this example, we are using the package numpy to create a series of x values ( Matplotlib hierarchy of figure elements, source . ::: -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, lines and text. Let's build up a simple line figure as an example. +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. -![matplotlib simple figure](images/mpl1.png){width="50%"} +![Matplotlib simple figure](images/mpl1.png){#fig-mpl1 width="50%"} -Note that - the behavior of `plt.show()` depends on how you run the script. If you are using Spyder and want `plt.show()` to work: - -1. Go to Tools. -2. Go to Preferences. -3. Select IPython console. -4. Go to Graphics tab. -5. In the Graphics backend section, select Automatic as the backend type. -6. Restart your kernel (go to Console > Restart kernel). - -Instead of a single plot, we can add different plots to the figure, for example two. Let's try to add another pot with the cosine values. +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. ```{python} # Create some data @@ -201,7 +193,7 @@ 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](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html) 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. +We can use the subplots method to create a figure and an array of two axes, one for each plot. Check the [subplots documentation](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html) 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. ```{python} import matplotlib.pyplot as plt @@ -212,7 +204,7 @@ f, axarr = plt.subplots(2, sharex=True) ``` ::: {.alert .alert-success} -> **Question 4**: `axarr` is an array. What are the elements of this array and how many elements does is exist out of? +> **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](https://matplotlib.org/stable/api/axes_api.html) to see what more you can tweak. @@ -223,11 +215,11 @@ 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 a two plots', fontsize=16) +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. +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. ```{python} # Axis label @@ -252,9 +244,7 @@ Finally! our plot is done for now. We have entered a lot of commands, stacked a plt.show() ``` -Note that if you paste the above code blocks cell by cell in a REPL, you will see the final figure as expected. In Jupyter Notebook, however, you may get nothing at the last 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 and 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. +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](images/mpl2.png){width="50%"} @@ -285,25 +275,25 @@ ax3.barh(x, y, color='y') plt.show() ``` -![matplotlib plot type examples](images/mpl3.png){width="50%"} +![Matplotlib plot type examples](images/mpl3.png){width="50%"} ::: {.alert .alert-success} -> **Question 5**: In the upper right subplot, why is there no point at x=3, y=8?. +> **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](https://matplotlib.org/stable/plot_types/index.html) and play around to find out more! -## Vizualizing spatial data +## Visualizing spatial data So, plotting sine and cosines is fun and all, but this is a course about **geo**scripting, 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`](https://github.com/SciTools/cartopy/blob/main/lib/cartopy/mpl/geoaxes.py#L354). 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](https://foundations.projectpythia.org/core/cartopy/cartopy.html) website, an excellent source for more information and tutorials about working with python! +For this tutorial, part of the material was taken from the [Project Pythia](https://foundations.projectpythia.org/core/cartopy/cartopy.html) 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. +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. ```{python} import matplotlib.pyplot as plt @@ -325,13 +315,13 @@ We don't see anything yet! That's because we have not put anything on the map. L 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](https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html) 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()` (don't believe it? [Check the source](https://github.com/SciTools/cartopy/blob/75939f9e81ac67838c52ce80230e4431badbeace/lib/cartopy/mpl/geoaxes.py#L609)! ) effectively. Have a look and play around with adding other features! Using the `linewidth` and`edgecolor` arguments we can style the map. Don't forget to `plt.show()` to see the map! +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](https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html) 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](https://github.com/SciTools/cartopy/blob/75939f9e81ac67838c52ce80230e4431badbeace/lib/cartopy/mpl/geoaxes.py#L609)! ). Have a look and play around with adding other features! Using the `linewidth` and`edgecolor` arguments we can style the map. Don't forget to `plt.show()` to see the map! ::: {.alert .alert-success} -> **Question 6**: Create a worldwide map with 3 different features, each styled differently. Also add the stockimage to the map. Use [the documentation](https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/geoaxes.html?highlight=stock#cartopy.mpl.geoaxes.GeoAxes.stock_img) to see how. +> **Question 5**: Create a worldwide map with 3 different features, each styled differently. Also add the stockimage to the map. Use [the documentation](https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/geoaxes.html?highlight=stock#cartopy.mpl.geoaxes.GeoAxes.stock_img) to see how. ::: -We have now step by step built up a map in a Pate Carree projection system. However, this projection has some major issues, have you seen [how big Antartica would be at the equator](https://www.thetruesize.com)?! Let's quickly create another map in another projection. In the [documentation](https://scitools.org.uk/cartopy/docs/latest/reference/projections.html) we can find a large list of projections that can be used. +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](https://www.thetruesize.com)?! Let's quickly create another map in another projection. In the [documentation](https://scitools.org.uk/cartopy/docs/latest/reference/projections.html) we can find a large list of projections that can be used. ```{python} fig = plt.figure(figsize=(11, 8.5)) @@ -344,47 +334,46 @@ 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 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. +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. ```{python} import geopandas as gpd gdf = gpd.read_file('https://raw.githubusercontent.com/GeoScripting-WUR/PythonProgramming/master/data/gadm41_NLD_2.json') -# Using the dutch coodinate reference system RDNew (epsg code 28992) -crs = ccrs.epsg(28992) +# Using the dutch coordinate reference system RDNew (epsg code 28992) +crs_rd = ccrs.epsg(28992) -# The data is in another projection as our plot, reprojection to RDnew +# The data is in another projection than our plot, reprojection to RDNew gdf = gdf.to_crs(28992) -# Initiate the plot, a little bigger then before +# Initiate the plot, a little bigger than before fig = plt.figure(figsize=(15, 15)) -ax = plt.subplot(1, 1, 1, projection=crs) +ax = plt.subplot(1, 1, 1, projection=crs_rd) ax.set_title('The municipalities of NL') -# Draw gridlines +# 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) +ax.set_extent((min_x, max_x, min_y, max_y), crs=crs_rd) -# ax.set_extent(gdf.total_bounds, crs=crs) would do this in one step, -# but the coordinates can be defined seperately as well in this order! +# 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, edgecolor = 'black', facecolor = 'None') +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](https://github.com/SciTools/cartopy/blob/main/lib/cartopy/mpl/geoaxes.py#L354) 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. +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](https://github.com/SciTools/cartopy/blob/main/lib/cartopy/mpl/geoaxes.py#L354) 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. ```{python} import contextily as cx @@ -403,11 +392,11 @@ The image tiles that are used as a basemap in this example are by default availa 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. ::: {.alert .alert-success} -> **Question 7**: By default, `contextily` uses OpenStreetMap tiles. Have a look at the [contextily documentation](https://contextily.readthedocs.io/en/latest/intro_guide.html#Providers) to find another tile provider, and use it instead. +> **Question 6**: By default, `contextily` uses OpenStreetMap tiles. Have a look at the [contextily documentation](https://contextily.readthedocs.io/en/latest/intro_guide.html#Providers) 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](https://github.com/rasterio/rasterio/blob/7c83de410b7b812e6552e4b76ce236bb6213c80d/rasterio/plot.py#L34) 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` `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. +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](https://github.com/rasterio/rasterio/blob/7c83de410b7b812e6552e4b76ce236bb6213c80d/rasterio/plot.py#L34) 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. ```{python} import requests @@ -421,25 +410,25 @@ 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, "data.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 = ccrs.epsg(32631) +crs_utm = ccrs.epsg(32631) # Initiate the area of interest fig = plt.figure(figsize=(15, 15)) -ax = plt.subplot(1, 1, 1, projection=crs) +ax = plt.subplot(1, 1, 1, projection=crs_utm) ax.set_title('The municipalities of NL') -# Read the GeoJson with GeoPandas. +# 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! +# 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 @@ -456,21 +445,95 @@ 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](https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.plot.html) 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. + +```{python} + +# 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](https://github.com/geopandas/geopandas/pull/3728). 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](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.colorbar.html) 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. + + +```{python} +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 short 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. +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 it's power. We made basic 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. +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! +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](https://docs.Python.org/3/contents.html) +- [Official Python tutorial](https://docs.python.org/3/contents.html) - [Python Style guide](https://www.python.org/dev/peps/pep-0008/) - [Python 3 Cheatsheet](https://ugoproto.github.io/ugo_py_doc/py_cs/) - [Overview Python package Cheatsheets](https://www.datacamp.com/community/data-science-cheatsheets?tag=python) diff --git a/pixi.toml b/pixi.toml index ef32426..e56c73c 100644 --- a/pixi.toml +++ b/pixi.toml @@ -10,7 +10,7 @@ version = "0.1.0" [dependencies] python = "*" cartopy = "*" -geopandas = "*" +geopandas = ">=1.1.4,<2" rasterio = "*" matplotlib = "*" numpy = ">=2.5.1,<3" @@ -24,7 +24,7 @@ channels = ["conda-forge"] [feature.python-programming.dependencies] cartopy = "*" python = "*" -geopandas = "*" +geopandas = ">=1.1.4,<2" rasterio = "*" matplotlib = "*" contextily = "*" diff --git a/stuff-i-removed-from-the-refresher.Rmd b/stuff-i-removed-from-the-refresher.Rmd deleted file mode 100644 index 53f1218..0000000 --- a/stuff-i-removed-from-the-refresher.Rmd +++ /dev/null @@ -1,378 +0,0 @@ -## Classes and objects - -Python is an *object-oriented programming language*. That is a programming paradigm that structures a code hierarchically with *classes* and *objects*. A class is a blueprint of functions and attributes to build an object, while an object is a self-contained component operationalizing the functions and attributes. - -For example, one could have a class `dog` with functions, `bark()` and `doginfo()`, and properties `breed` and `age`. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -class Dog: - def __init__(self, breed, age): - self.breed = breed - self.age = age - - def bark(self): - print("bark bark!") - - def doginfo(self): - print("This " + self.breed + " is " + str(self.age) + " year(s) old.") -``` - -One can then create an *instance* of this class, i.e. an object, that represents Ozzy, a Maltese of two years old. This object can apply the function `eat()` and `run()`. We can call the attributes and functions with the dot notation: - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -ozzy = Dog("Maltese", 2) -print(ozzy.breed) -print(ozzy.age) -ozzy.bark() -ozzy.doginfo() -``` - -During Geo-scripting, you do not need to define your own classes. However, sometimes existing classes of modules will be used, such as a `DataFrame` in Pandas or a `GeoDataFrame` in GeoPandas. It is therefore important to understand the concepts of classes and objects generated from them. - -```{block, type="alert alert-success"} -> **Question 1**: Suppose that we want to build a script in Python of an orchestra playing a song. What would be the class and what would be the objects? -``` - -## Data types and variables - -Values belong to a data type. The most important built-in (base) types are: -* integer -* floating point -* Boolean -* several compound data types (see next section), like string - -Python is a strongly typed language. That is, you cannot perform operations inappropriate to the type of the value. For example, attempting to add integers to strings will fail. For this reason, understanding and being aware of value types in Python is crucial. - -Values can be cast to other types, e.g.: - -```{Python,engine.path='/usr/bin/python3'} -print(int(10.6)) -``` - -```{block, type="alert alert-success"} -> **Question 2**: What is the difference between 10 and 10.0 in Python? -``` - -Often, we do not use values directly in a script. Instead, we use access them through variables. A variable is a way to reference to a known or unknown value. In Python, a value is assigned to a variable with the `=`-sign, e.g.: - -```{Python,engine.path='/usr/bin/python3'} -building = 'Gaia' -buildingnumber = 101 -print(building + ' is in Wageningen') -``` - -As opposed to other programming languages, in Python, the data type does not need to be explicitly defined when creating a variable. Python derives the data type of the variable from the value assigned to it. In the code above, the variable `building` is thus a string, and the variable `buildingnumber` is an integer, without having defined this explicitly. - -Variable names should: 1) start with a lowercase letter, and 2) not contain spaces. Furthermore, it is advisable to use meaningful variable names, both for yourself and for (future) users of your code. - - -## Compound data types - -*Compound data types* or container types are Python data types that have in common that they can be broken down into smaller *elements*. The most commonly known compound data types are: -- string -- list -- dictionary -- tuple -- NumPy array - -For most compound data types, the elements can be accessed by providing the index of the element with the square-bracket operator, `[ ]`. This works, for example, for items in lists and letters in strings. Notice that indexing in Python, as opposed to other languages like R, **STARTS AT 0!** Positive indexes access elements from the front, and negative indexes from the back. Multiple elements can be selected with a colon `from:to`, where the first is inclusive and the last exclusive. Here is an example for accessing items in a list, and in a string: - -```{Python,engine.path='/usr/bin/python3'} -# List -campus = ['Gaia', 'Lumen', 'Radix', 'Forum'] -print(campus[3]) # Forum -print(campus[-1]) # the last item of the list -print(campus[0:3]) # the first 3 items (index 0, 1, and 2), i.e. 'slicing' -our_building = campus[0] # GAIA - -# String -print(our_building[0]) # G -print(our_building[0:2]) # GA -``` - -```{block, type="alert alert-success"} -> **Question 3:** What building is `campus[-2]`? Test it. -``` - -Lists can be nested. That means that a list is an item in another list. Accessing items in nested lists works the same as in regular lists, e.g,: - -```{Python,engine.path='/usr/bin/python3'} -# Nested list -samples = [["x", "y", "z"], [12, 32, 7], [12, 40, 7]] - -# Access the x -header = samples[0] -first_item = header[0] - -# Or at once -first_item = samples[0][0] -``` - - -```{block, type="alert alert-success"} -> **Question 4:** How can we access the value 40 above? Test it. -``` - -A *dictionary* is set of key-value pairs. The key can be used to access the value. A dictionary is different from the other compound data types as it is not ordered, i.e. the order of the elements is undetermined. As such, the square-bracket operator, `[ ]` cannot be used. - -```{Python,engine.path='/usr/bin/python3'} -# Dictionary -campus_dictionary = {101: 'Gaia', - 100: 'Lumen', - 107: 'Radix', - 102: 'Forum', - 104: 'Atlas'} - -# Access dictionary value using key -print(campus_dictionary[102]) -``` - -In contrast to a list, a *NumPy array* is a *homogeneous multidimensional array* from the NumPy package. Homogeneous refers to the fact that all data in one array have to be of the same type. NumPy’s array class is called `ndarray` (n-dimensional). In the array, the dimensions are called *axes*, and the number of axes is the *rank*. - -
-Anatomy of a numpy array -
Anatomy of a NumPy array, source https://valecs.gitlab.io/resources/numpy.pdf.
-
- -Below, we show how to create a NumPy array from a set of lists. We also show how to create *standard* arrays, filled with zeroes, ones, or random number, often to be updated with actual data later in your script. - -```{Python,engine.path='/usr/bin/python3'} -import numpy as np - -# Create array from list -a = np.array([[1, 3, 4], [2, 7, 6]]) -print('a is', a) - -# Create standard arrays. -print(np.zeros((3, 2))) -print(np.ones((2, 3), int)) -print(np.ones((2, 3), int) * 5) -print(np.empty((2, 2))) -``` - -Just like with other compound data types, the elements can be accessed with the square-bracket operator, `[ ]`. But in NumPy arrays, there is one index per dimension (per axis), separated by commas. Multiple elements per dimension can be selected with a colon (`from:to`, just like in a list). Leaving out the index before or after the colon, selects all elements from the beginning or to the end, respectively, see the examples in the figure below. - -How to slice NumPy arrays - -Finally, *structured NumPy arrays* can have a different data type per column, i.e. per attribute. It is thereby similar to a `DataFrame` in R (we'll look at another similar data type tomorrow). In structured NumPy arrays, you can access or edit attribute values either by dimension or by attribute name. - -```{Python,engine.path='/usr/bin/python3'} -# Create create structured NumPy array filled with zeroes -data = np.zeros(4, dtype={'names': ('name', 'age', 'weight'), - 'formats': ('U10', 'i4', 'f8')}) -print(data.dtype) - -# Now we can fill this structured array with data (lists) of the correct type -name = ['Alice', 'Bob', 'Cathy', 'Doug'] -age = [25, 45, 37, 19] -weight = [55.0, 85.5, 68.0, 61.5] -data['name'] = name -data['age'] = age -data['weight'] = weight - -# Inspect the result -print(data) -``` - -## Functions - -A function is a section of code with a common purpose. Functions can be useful for 1) making the main part of the code concise, 2) debugging, as functions can be tested and edited separately from the main code, and 3) reuse of code in the same or other programs. - -Functions may (not always) accept *argument(s)*, perform an *action*, and may (not always) *return value(s)*. For example, the built-in function 'int()' we have seen before performs the action of converting a value to integer type: - -```{Python,engine.path='/usr/bin/python3'} -# Return value into a variable, the function name, its (single) argument -myint = int(10.6) -``` - -You can also create functions yourself. The syntax for creating a function is: - -``` -def functionname(arg1, arg2, argn): - ... (expression1) ... - ... (expression2) ... - ... (expressionn) ... - return returnvar1, returnvar2, returnvarn - -``` - -For example, we could create a function to compute the area of a rectangle based on its width and length, provided as arguments, and return the output: - -```{Python,engine.path='/usr/bin/python3'} -def calculate_rectangle_area(width, length): - rectangle_area = width * length - return rectangle_area - - -print(calculate_rectangle_area(4, 3)) - -# Or alternatively; -output_number = calculate_rectangle_area(width=4, length=3) -print(output_number) -``` - -Functions or classes can be made more informative by *docstrings*: - -```{Python,engine.path='/usr/bin/python3'} -def calculate_rectangle_area(width, length): - """Computes the area of a rectangle by multiplying width and length. - :width: width is the width of the rectangle - :length: length is the length of the rectangle - :returns: width * length - """ - rectangle_area = width * length - return rectangle_area - - -print(calculate_rectangle_area(4, 3)) -``` - -Variables created in functions are *local*. That means that a variable created in a function does not exist outside of the function. So, here, `rectangle_area` cannot be accessed outside of the function definition (after the indented part of the code). - -```{block, type="alert alert-success"} -> **Question 5**: Test typing `print(rectangle_area)` at the end of the code block defined above. What happens? Why is that? -``` - -## Modules and packages - -Python packages and modules are collections of classes and functions. Any Python file is a _module_, its name being the file's base name without the `.py` extension. A _package_ is a directory of Python modules containing an additional `__init__.py` file, to distinguish a package from a directory that just happens to contain a bunch of Python scripts. Packages can be nested to any depth, provided that the corresponding directories contain their own `__init__.py` file. - -Basically there are four ways to load a module and/or a package: -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -import math -print(math.pi) - -from math import pi -print(pi) - -from math import pi as ip -print(ip) - -import numpy as np -print(np.pi) -``` - -```{block, type="alert alert-success"} -> **Question 6**: What are the differences between these four ways to import modules and/or packages? -``` - -Often-used internal packages/modules: - -- `os`: Operating system features -- `sys`: System specific configuration -- `math`: Mathametical functions, operators and constants -- `datetime`: Date/Time functionality - -## Conditionals - -*Comparison operators* compare two values or, more commonly, variables. The operands (x and y below) are often numerical, e.g. floating point or integer. The result of comparison operators is a 0 (False) or 1 (True), of type Boolean. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -x == y # TRUE if x is equal to y -x != y # TRUE if x is not equal to y -x > y # TRUE if x is greater than y -x < y # TRUE if x is less than y -x >= y # TRUE if x is greater than or equal to y -x <= y # TRUE if x is less than or equal to y -``` - -*Logical operators* evaluate the logical relation between two values or variables. The operands (x and y below) are in most cases Boolean. The result of logical operators is also Boolean. - -```{Python,engine.path='/usr/bin/python3', eval=FALSE} -x and y # TRUE if both x and y are TRUE -x or y # TRUE if x or y are TRUE -not x # TRUE if x is FALSE -``` - -A *conditional statement* checks whether a condition is fulfilled and only if it does, it executes a block of code. The syntax of a conditional statement is: - - -``` -if condition: - ... (expression1) ... - ... (expression2) ... - ... (expressionn) ... -else: - ... (alternative_expression1) ... - ... (alternative_expression2) ... - ... (alternative_expression3) ... - -``` - -The condition should be a Boolean variable, or an expression resulting in a Boolean. A condition is typically built up using comparison operators and/or logical operators. The expressions after the `if`-statement are executed when the condition is met, while the conditions in the `else`-statement (not mandatory to include) are executed when the condition is not met. Multiple conditions can be checked consecutively with one or more `elif`-statements in-between the `if` and `else`. For example: - -```{Python,engine.path='/usr/bin/python3'} -x = 3 -if x == 1: - print("it is one") -elif x == 2: - print("it is two") -elif x == 3: - print("it is three") -else: - print("above 3") -``` - - -```{block, type="alert alert-success"} -> **Question 7**: Think of a conditional statement that uses logical operators in the condition instead of comparison operators. - -``` - -## Loops - -Loops are an essential construct in Python, more than in R. Two types of loops exist: -- the `for`-loop, used when it can be known beforehand how many iterations are required -- the `while`-loop, used when it cannot be known beforehand how many iterations are required - -The `for`-loop is often the preferred construct, as it is computationally more efficient and less error-prone than the `while`-loop; no chance to end up in an endless loop or to accidentally skip an element (iteration). The syntax of a `for`-loop is: - - -``` -for element in compound: - ... (expression1) ... - ... (expression2) ... - ... (expressionn) ... - -``` - -Often `i`, standing for iterator, is used as the `element`, but it can be any variable name. Key is that the current element is put into that variable in each iteration, so the `element` is what one should use in the expressions in the `for`-loop. The `for`-loop works on any compound data type: - -```{Python,engine.path='/usr/bin/python3'} -# Looping over a list -campus = ['Gaia', 'Lumen', 'Radix', 'Forum'] -for i in campus: - print(i) - -# Looping over a string -our_building = campus[0] # GAIA -for letter in our_building: - print(letter) -``` - -The `while`-loop works on a condition instead of on a compound data type. The syntax of a `while`-loop is: - - -``` -while condition: - ... (expression1) ... - ... (expression2) ... - ... (expressionn) ... - -``` - -In a `while`-loop, when an element is used, it is not automatically altered in each iteration; the alteration needs to be explicitly programmed. - -```{Python,engine.path='/usr/bin/python3'} -n = 0 -while n < 20: - print(n) - n = n + 1 -``` - - -```{block, type="alert alert-success"} -> **Question 8**: What is the last number printed? Now print n after the loop (not in the indented part); what does it print? -``` - -You may remember that variables in a function are local. This is different for loops: Variables in a loop are NOT local, and thus affect what happens outside of the loop, as demonstrated by the last question. Be aware of this. \ No newline at end of file diff --git a/styles.css b/styles.css index 6dd1c79..35954e6 100644 --- a/styles.css +++ b/styles.css @@ -34,3 +34,10 @@ td code:not(.sourceCode) { vertical-align: middle; } +/* Quarto only centers the image

/

inside a centered figure, not + the
itself, so a resized image's caption is left-aligned + below it instead of centered under the image. */ +.quarto-figure-center figcaption { + text-align: center; +} +