Multispectral Satellite Imagery: Processing Using Python

Landsat-8

Landsat-8 is a multispectral imaging satellite launched by NASA. The satellite has a 16-day repeat cycle. Other specifications are given below:

 

Landsat-8 images can be downloaded from https://earthexplorer.usgs.gov/.

After downloading data, we can start writing code to perform operations. Here we will learn how to read raster file, layer stack it, to read shapefile of the given area and mask the raster as per shapefile extents.

The first step is to import the required packages. Rasterio provides functionality for reading, writing and performing operations on raster data formats. GeoPandas is used to perform spatial operations on geometric data types.

import rasterio
from rasterio import mask
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd

Then assign the path of all 11 bands.

b1 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B1.TIF' 
b2 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B2.TIF' 
b3 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B3.TIF' 
b4 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B4.TIF' 
b5 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B5.TIF' 
b6 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B6.TIF' 
b7 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B7.TIF' 
b8 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B8.TIF' 
b9 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B9.TIF' 
b10 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B10.TIF' 
b11 = 'LC08_L1TP_148047_20190325_20190403_01_T1_B11.TIF'

It is not convenient to read different bands and process it seperately. Instead, we can stack all bands in a single raster file for further processing. We can use the following code for layer stack operation. The stacked raster file is saved on the given path.

file_list = [b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11]
stacked_path = 'stacked_raster.tif' #Path of new stacked file

# Read metadata of first file
with rasterio.open(file_list[0]) as src0:
    meta = src0.meta

# Update meta to reflect the number of layers
meta.update(count = len(file_list))

# Read each layer and write it to stack
with rasterio.open(stacked_path, 'w', **meta) as dst:
    for id, layer in enumerate(file_list, start=1):
        with rasterio.open(layer) as src1:
            dst.write_band(id, src1.read(1))

We can read stacked raster using rasterio built-in function.

raster = rasterio.open(stacked_path)

If one needs to check information about projection and geo transform of input raster, this can be done as follows:

raster.crs

Out:

CRS.from_epsg(32643)

raster.transform

Out:

Affine(30.0, 0.0, 110385.0,
       0.0, -30.0, 2196615.0)

Now for reading raster as an array:

raster_arr = raster.read()
print(raster_arr.shape)

Out:

(11, 7841, 7681)

Create FCC (False Color Composite) for visualization:

If we represent a multispectral image using wavelength regions other than visible red, green and blue, it is called false-color composite (FCC). By using FCC, data interpretability can be enhanced as spectral separation is increased. Hence we create FCC for visualizing the image.

raster_fcc=np.dstack((raster_arr[4,:,:],raster_arr[3,:,:],raster_arr[2,:,:]))

We need to scale down values to display using matplotlib,  a python library.

fig = plt.figure(figsize = (15,15))
plt.imshow(raster_fcc/32768)
plt.title('FCC543')

Input shapefile

The administrative area shapefile of India can be found at https://www.diva-gis.org/gdata. We use geopandas built-in function for reading shapefile. We will clip raster as per shapefile of the Mumbai city of Maharashtra, India.

ind_gdf = gpd.read_file('IND_adm/IND_adm2.shp')
ind_gdf.head()

state_gdf = ind_gdf[ind_gdf['NAME_1'] == 'Maharashtra']
state_gdf.head()

city_gdf = state_gdf[state_gdf['NAME_2'] == 'Greater Bombay']
city_gdf

The next step is to clip raster using a shapefile of the given area. Clipping is a process of masking out regions of a raster that fall outside the boundary of polygons defined in the shapefile.

For performing this operation, raster and shapefile both should be in the same CRS. Here we use ‘32643’ as EPSG of UTM zone.  This can be changed if the geographic area is changed.

city_gdf = city_gdf.to_crs(epsg=32643)

For clipping, we will use mask method which takes raster file and vector file datasets as inputs

raster_path = stacked_path
with rasterio.open(raster_path) as src:
    out_image, out_transform = rasterio.mask.mask(src, city_gdf.geometry, crop=True)
    out_meta = src.meta

Now create FCC array for visualization.

city_fcc = np.dstack((out_image[4,:,:], out_image[3,:,:], out_image[2,:,:]))

Again, the values have been scaled down to display using matplotlib.

fig = plt.figure(figsize = (15,15))
plt.imshow(city_fcc/32768)

 

Finally,  we get the raster file of the given city clipped using the given shapefile.

References

1. Jensen, John R. Remote sensing of the environment: An earth resource perspective 2/e. Pearson Education India, 2009.
2. https://www.usgs.gov/media/images/landsat-8-band-designations
3. https://earthexplorer.usgs.gov/
4. https://www.diva-gis.org/gdata
5. https://rasterio.readthedocs.io/en/latest/topics/masking-by-shapefile.html
6. https://eos.com/false-color/

For any queries related to this blog, you can reach me at LinkedIn