This is the second post in a series that will teach non-programmers how to use Python to handle and analyze geospatial data.
Introduction
This is the second story of the series Python for Geosciences — working with satellite image data. In the first post of this series (here) we set up the environment to run Python code from a Jupyter Notebook and learned how to open a GeoTiff image by using the rasterio package. Today we will learn the basics of matrix manipulation. This will allow us to combine matrices to create data cubes, perform raster calculations and create spectral indices.
As our series is based in step by step examples, our goal is to calculate the Modifies Normalized Difference Water Index (MNDWI) and select just the pixels with a high probability of being water. At the end, we will display the water mask and the RGB image side by side for comparisons. I will write the code throughout the text as snippets, but at the end I will provide the full notebook.
Step 1- Open the image (several bands)
Similarly to the previous post, we will open the downloaded Landsat image from Manaus/Brazil region. However, instead of opening one file we need to open several bands that are needed to calculate the index and the final RGB image. Probably the more straightforward way (and that you will be tempted to) to open the various images would be to repeat each line of the previous post and store each band as a different variable, as bellow:
b2 = rasterio.open(#image_name_here)
b3 = rasterio.open(#image_name_here)
b4 = rasterio.open(#image_name_here)However, in order to write a cleaner and reusable code, it is better to write a specific function for it and “organize” all the image bands into a single object. Additionally, we will learn how to deal with dictionaries and use the dictionary comprehension syntax.
Python Dictionary
A Python dictionary is an object that can store pairs of keys and values. Keys can be of any “internal” Python’s data types (i.e. string, numbers or tuples) and values can be of any arbitrary type/object (even an image). Python dictionaries are declared with the following syntax:
variable = {key: value, key: value2, …}
Note: Until Python 3.6, the order of the items within a dictionary were not guaranteed to remain the same. We said that the Python’s dictionary was unordered. That has changed in Python 3.7, where the key:value pairs remains in the same order they were initially declared.
So, the main idea will be to have the bands of our image declared into a single object like this:
img = {'B1': band1, 'B2': band2, 'B2': band3, ...}Note that the keys ‘B1’, ‘B2’ and ‘B3’, are defined arbitrarily so we can access our data later, like so: img[‘B1’] or img[‘B2’]. Moreover, band1, band2 and band3 could be any object we want. For example, we could store the paths for the files in the filesystem or the raster image already loaded in memory. It is really a project’s design decision and depends on our necessity. For example, we could write:
img = {
'B1': 'D:/landsat/1026_20201106_02_T1_SR_B1.TIF',
'B2': 'D:/landsat/1026_20201106_02_T1_SR_B2.TIF',
...}Or,
band1 = rasterio.open('D:/landsat/...02_T1_SR_B1.TIF').read(1)
band2 = rasterio.open('D:/landsat/...02_T1_SR_B2.TIF').read(1)
img = {'B1': band1, 'B2': band2}In the first example, our image stores not the bands themselves, but the path of each GeoTiff file whereas in the second example our image stores the actual bands values. In our case, we will get the second option and define a simple function to open the bandas we want.
Defining a function
In the first post we used functions already defined in the packages we imported such as open and read from rasterio and the imshow from the matplotlib. Now it is type to write out own function.
A function is a piece of code that receives some arguments as inputs and returns an output (after some computations) and can be reused many times in the code, or can be reused by other codes. It is very simple to define a new function in Python. Put simply, it has the following syntax:
def name_of_the_function(arg1, arg2, arg3=1):
v = arg1+arg2/arg3
return vIn this syntax example, arg1 and arg2 are mandatory arguments and should be passed to the function, otherwise it will raise an error, whereas arg3 is an optional argument that assumes the default value 1 if it is omitted from the function call. To call this function, one should only write its name and pass in parentheses the necessary args.
When calling a function, there exists two different ways to pass the arguments to the function. We have positional arguments, where they are passed without a name and considered in the order they are passed and keyword arguments, that comes with the name of the argument preceding its value. For example:
result = name_of_the_function(10, arg2=3)In this example, arg1 was passed as a positional argument and arg2 was passed as a keyword argument. More on this topic can be found in https://www.programiz.com/python-programming/function-argument or other online material. Be careful, though, as there are some confusion explanation that mixes the concepts. Positional and keyword arguments are concepts related just to function calls, and not function definitions. This StackOverflow thread tries to clarify it (here). Now that we know how to define a new function, we are good to create our loading function.
The Load Image function
The following code creates our function to load a Landsat image.
This functions has some points that are worth mentioning. The first one, and probably the scariest for beginners, is the next(path.glob(f’*_SR_{band}.tif’)).
To understand this line of code, we have to split it up in smaller pieces. First we will take a look at the string. The letter f before the string indicates that it is a Python’s f-string. This string can receive variables inside curly brackets and outputs everything as a single string. In our case we want strings like “*_SR_B2_.tif” or “*_SR_B3_.tif” to serve as pattern strings to find the files we want, where * indicates everything. In our case, everything that ends with “_SR_B2.tif” or “_SR_B3.tif” will be matched.
Path.glob is the function to search for files that match a specific pattern (glob is the short for global and if you are curious about it, you can check it here. The result of a glob function is usually a list of files that matches the pattern. In our case, we just want the first file (it will hopefully exist just one) so we will use the next function.
Note: The correct explanation for the next(path.glob) passes through the understanding of generators in Python. It is not complex (if corrected explained), but we will keep things simple for now.
Step 2- Accessing image values
With the image loaded in memory and organized in the dictionary, we can access the bands by calling the image variable followed by the keys inside square brackets, like so: img['B1'] or img['B2'], etc. The type function displays the data type of a specific variable. Note that our bands are arrays (a data type defined in the NumPy package).
type(img), type(img['B2'])
code output:
(dict, numpy.ndarray)In the notebook I’ve added a cell to display the 3 loaded bands side by side (Figure 1).
Display image as RGB
The first thing we will do with these data is to display it as a natural RGB colors. From Landsat documentation we know that B4=Red, B3=Green and B2=Blue. Additionally, if we look at the imshow documentation (the function we’ve used to display the image), we will see that the X argument can be of shape (N, M, 3), for an image with RGB values (0–1 float or 0–255 int). Our bands are already loaded with (N, M) matrices, so we need to combine them in the correct order to create a “Cube” and show it correctly (Figure 2).
The NumPy command that will use to combine the bands is the stack. One important remark has to be done concerning the axis parameter. The default behavior of the stack is to index the joined matrices in the first axis (considered axis=0). That would make a final cube with shape (3, N, M). But the imshow needs the RGB channel indexed in the last axis. In this case we will use the parameter axis=-1.
Display RGB Function
Just as before, we will write a simple function to display an image as RGB considering the bands names that we pass in. This way, we will be able to create false color images using other bands.
Within the function, I added a line to divide the values by the maximum value, because the imshow requires the values to be in the range [0..1]. As this image has clouds that are over bright compared to the land, I added an argument to increase the alpha (brightness) just for visualization purposes. The internal values of the bands are not changed at all. Figure 3 shows the same image shown in natural colors and in the false color composite of the last line of code.
Step 3- Create a Normalized Difference index
Now that we know how to open and display the bands of the Landsat image, we can compute the spectral index we want. For that, we will perform element wise operations using different bands. Element wise operations can be performed in matrix by using built-in math operators (+, -, * or /) or other NumPy functions, such as square, sqrt, and others. For example, adding the values of two bands (considering our image as dictionary) is as simple as:
addition = img['B1'] + img['B2']In this example, the result is another matrix with the element-wise addition of the two bands. We can then write a universal Normalized Difference (ND) function that can perform any ND index. A ND index is the difference divided by the sum (formula bellow) and its values can range from [-1..+1]:
nd = (Band1 - Band2) / (Band1 + Band2)The following code shows the function to calculate a ND index. The two first lines of code are there to ensure we don’t have division by 0, that is something can occurs with both values happens to be 0. In this case, when both are 0 at the same time we say it is NaN (Not a Number). Note that the equals (==) and the AND (&) operator are also element wise, so they return a matrix of Boolean (Trues or Falses). We will explore more of this topic when performing matrix slicing.
The output of the code can be seen in Figure 4. Note that the NDVI index gives a better response (higher values = greener because of the colormap we applied) over the vegetation areas and the MNDWI has a higher response (higher values = bluer because of the ‘Blues’ colormap) in the rivers.
Figure 4: (a) NDVI index displayed in a “greens” colormap, where higher values are represented as greener pixels and (b) MNDWI index displayed in a “blues” colormap where higher values are represented as bluer pixels. Image by author.
Step 4- Masking by index values
Now that we have both indices stored in memory, suppose we would like to select just the water portion of the image, or just the vegetation portion. For that, we have to understand array slicing and masking. As it is a long topic that must be very well mastered by those working with geospatial images, we will not cover everything on this issue. We will see just the basics today and next post we will do an in-depth analysis on it, so, don’t be in panic if you don’t understand all the details today. And don’t hesitate to leave in the comments any question or demand that we could discuss in the next story.
When we perform an element-wise boolean operation such as equal (==), different (!=), greater than (>), less than (<), AND (&) or OR (|) to keep up with the most commons, the result is a boolean array. So, for example, if we type:
img['B3'] > img['B5']The result will be a boolean array with the same shape of the bands (remember that this is an element-wise operation, so the shape doesn’t change) with TRUE when B3>B5 and false otherwise. So, if we want to select the portion of the image that has “water”, we select arbitrarily a threshold for the MNDWI index, above which we will consider water. Water identification through satellite images is much more complex than that, and depends on many other factors. You can check my story Water Detection in High Resolution Satellite Images using the waterdetect python package to know more about it.
In this tutorial, for the sake of simplicity, we will consider water pixels those that have MNDWI value greater than 0.0. If we type water_mask = mndwi> -0.1 we will receive a Boolean array with TRUE on the “water” pixels and FALSE otherwise. We can then display this mask as any other array using the imshow function.
water_mask = mndwi > 0.0
water_mask
result:
array([[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
...,
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False]])
plt.figure(figsize=(7,7))
plt.imshow(water_mask)
Notebook
The full code to run these examples is provided in the Notebook bellow.
Other Posts
Link for other posts of Python for Geosciences series:
Part 1- Working with Satellite Images
Part 3- Spectral Analysis
Part 4- Raster bit masks explained
Part 5- Raster Merging, Clipping and Reprojection with Rasterio
Part 6- Scatter plots and PDF Reports
Conclusion
Today we learned the basics of matrix operations and we gave the first steps on manipulating dimensions, stacking layers to create composites. If you feel comfortable to go deeper, you can try to create other indices and use them to display false color with the indices.
Other option, that we will learn in the next story is to separate the actual values of water pixels or to overlay the mask on a RGB image, like in the examples of Figure 6.
So, stay tuned and see you in the next story!