Vincent GODARD - V1 - 29/01/2024
Freely inspired by :
Sources :
Landsat Thematic Mapper (TM) data from September 10, 1987, West of Worcester, Massachusetts (from the https://clarklabs.org/download/)
scikit-image https://scikit-image.org/
matplotlib https://matplotlib.org/
Compressed file to download => here.
Upload the following file from your storage area to the data directory (TD13/data) under this notebook is located (.ipynb):
○ "how87tm432.jpg" ;
○ "how87tm321_PS.jpg"
This noteboock and the "data" directory will be at the same level in the tree, both under TD13.
Some libraries are already preinstalled (e.g. Matplotlib). The Skimage library is not, and must be installed.
You then need to load them (Import function).
# To install Skimage, delete the # (before the "!").
# !pip install scikit-image
## Loading libraries
import skimage # Image processing library
import skimage.io as io # Utility for reading and writing images in various formats (shortened to io).
import skimage.metrics # Utility allowing measurements between images.
import matplotlib.pyplot as plt # Library for plotting and visualizing data in graphical form (shortened to plt).
import numpy as np # Library intended to manipulate multidimensional matrices or tables.
## To find out the version of skimage
skimage.__version__ # warning to "." after skiimage
Importing the skimage version may generate an error message, such as incompatibility with numpy, for example. One of the libraries is too old for the other. An update is therefore required.
## SciPy is a project that federates Python libraries for scientific use.
## Scipy uses arrays and matrices from the NumPy module (https://fr.wikipedia.org/wiki/SciPy).
## Update scipy, delete the # (before the "!") if necessary.
#!pip install --upgrade scipy
## Next, import it
#import scipy
## Then check that this new version is compatible
#scipy.__version__
You may find it useful to take a look at the https://scikit-learn.org/stable/ or https://scikit-image.org/ libraries, some of which have their own wikis in French^^! : https://fr.wikipedia.org/wiki/Scikit-learn or https://fr.wikipedia.org/wiki/Scikit-image or even https://fr.wikipedia.org/wiki/NumPy.
Reading an image file centered on the Howe Hill region, West of Worcester in Massachusetts (USA).
This is a natural color composite of the Landsat Thematic Mapper (TM) bands TM3 (red band), TM2 (green band) and TM1 (blue band) from September 10, 1987.
TM3 (red band) will be red-coded, TM2 (green band) will be green-coded and TM1 (blue band) will be blue-coded.
## Reading with the imread() function before displaying
HOW87TM321 = io.imread('data/how87tm321_PS.jpg')
The dynamic range of the image has been improved beforehand, so that it doesn't look too dull when first displayed, in this case with Photoshop (hence the _PS). We'll see later how to do this in Python.
## Displaying the number of rows, columns and number of channels (if applicable)
## Using the .shape function of the .io module
print(HOW87TM321.shape)
## If you want to dress up the output (text between '', then "," before the .shape function on the HOW87TM123 variable)
print('(number of rows, number of columns, number of bands):', HOW87TM321.shape)
We may also want to know the depth of our variable (its coding) by its type.
○ Classes are a way of bringing together data and functionality: https://docs.python.org/fr/3/tutorial/classes.html
○ The type of a Python object determines what kind of object it is: https://docs.python.org/3/glossary.html#term-type## Displaying information
print( 'classe :', type(HOW87TM321) )
print( 'type :', HOW87TM321.dtype )
The class is a numpy.ndarray vector placed in the HOW87TM321 image variable.
As it has 3 channels (cf .shape), it's a 3-dimensional array, so it's a 3D array of type NumPy.
The type is an integer (unsigned integer 8bits) encoded on one byte (256 levels of gray, for example).
So, our HOW87TM321 image is a vector of class numpy.ndarray containing values of type uint8 (single-byte integer) in 3 dimensions, whose form is: nb rows, nb columns, nb channels coded in Red, Green, Blue.
## Display the HOW87TM321 color composite image, using the imshow() function in the io module:
io.imshow(HOW87TM321)
The imshow() function in the io module calls matplotlib!
Display by iteration (Loop) using the number of channels (see Contents in 2.2).
## Decomposition into 3 RGB channels displayed in standard colors
## With a loop that traverses the range of channels "c" a specified number of times (3), starting from the first channel.
## Watch out for the ":" at the end of the for loop.
for c in range(3):
plt.imshow(HOW87TM321[:,:,c]) # [:,:,c] would allow under-window cutting (cf. SEMDOC/TD12)
plt.show() # displays images queued by imshow() the line before.
## To see examples of 'for' loops => https://waytolearnx.com/2019/05/boucle-for-en-python.html?utm_content=cmp-true
To the eye, there's no obvious difference between the channels! This will be more pronounced with PIR, band 4, for the false color image (see chap. 4. below)!
## Decomposition into 3 RGB channels displayed in gray and centered on a sub-zone
for c in range(3):
plt.imshow(HOW87TM321[400:480,200:300,c], 'gray') # thumbnails centered on the airport (approx. coordinates read in 3.1)
plt.show()
To the eye, the difference in grayscale is still not obvious!
Let's record these 3 channels and compare them.
If you're only interested in one channel, there's no need to loop.
Remember: channels are counted from 0!
0 = channel 3 (red band), 1 = channel 2 (green band) and 2 = channel 1 (blue band)
## Select the channel to be recorded, here 3.
io.imsave('data/HOW87TM3.png', HOW87TM321[:,:,0])
If you want to record all three channels, it's best to use a loop.
## Creation of a band number equivalent to the range increment
## We create a tuple, a list, where two pieces of information are linked and not permutable.
img_index_tuple_list = [(0,'TM3'),(1,'TM2'),(2,'TM1')]
## Successive selection of channels to be recorded with band number assignment
for c in img_index_tuple_list:
out_img_name = f'data/HOW87{c[1]}.png' # We tell it where to start [1], in the tuple
io.imsave(out_img_name, HOW87TM321[:,:,c[0]]) # incremental recording step
print(out_img_name) # to see below what the file names will be in /data
There doesn't seem to be much "visual" difference between the three .png files created. It would be better to check their differences by comparing their histograms, for example. This will be the subject of the next tutorial. In the meantime, you can test what you've just done on a new, higher-contrast color composition. A "false color" composition.
As we saw in TD 12, we can compare channels two by two using the modulus of mean square error (MSE).
Loading images into memory.
## Reading images with the imread() function before displaying them
HOW87TM3 = io.imread('data/how87tm3.png')
HOW87TM2 = io.imread('data/how87tm2.png')
HOW87TM1 = io.imread('data/how87tm1.png')
We can make sure that the images are the same size. This is a prerequisite for the two-by-two comparisons in the next paragraph (calculating MSE).
## Display number of rows, columns and channels (if applicable)
## Use the .shape function in the .io module
print(HOW87TM3.shape)
print(HOW87TM2.shape)
print(HOW87TM1.shape)
## Calculation of the mean square error (MSE) between the two channels TM3 and TM2
mse = skimage.metrics.mean_squared_error(HOW87TM3, HOW87TM2) # calculation of the gap between the two channels with skimage.metrics via MSE
MSE = f'{mse:.2f}' # MSE value to two decimal places, magnitude without unit (apparent luminances or numerical counts)
## MSE display
print('MSE (mean square error) value between TM3 and TM2 channels:', MSE)
Although the difference in radiometries is not obvious to the eye, the MSE value is unmistakable. The channels are different.
This can still be tested between TM3 and TM1.
## Calculation of the mean square error (MSE) between the two channels TM3 and TM1
mse = skimage.metrics.mean_squared_error(HOW87TM3, HOW87TM1) # calculation of the gap between the two channels with skimage.metrics via MSE
MSE = f'{mse:.2f}' # MSE value to two decimal places, magnitude without unit (apparent luminances or numerical counts)
## MSE display
print('MSE (mean square error) value between TM3 and TM1 channels:', MSE)
The TM3 and TM1 channels are closer than the TM3 and TM2 channels because their MSE is lower.
Can you calculate the MSE between channels TM1 and TM2?
Can you display another CC, HOW87TM432? Break it down into three bands?
The differences will be more visible on this "false color" composite!