Earth observation satellites collect data using two methods, passive and active.
Passive sensors use the sun as an illumination source. Sun’s energy is either reflected, absorbed or re-emitted. The electromagnetic radiation which is reflected or emitted from the earth’s surface is recorded by passive sensors. Examples include multispectral and hyperspectral sensors.
On the other hand, active sensors use their own illumination source and record the backscattered energy. Active sensors are used in Microwave and LiDAR remote sensing.
Landsat-8 is a multispectral imaging satellite launched by NASA. The satellite has a 16-day repeat cycle. Other specifications are given below:
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) 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))
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:
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)
(11, 7841, 7681)
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')
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.
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
city_fcc = np.dstack((out_image[4,:,:], out_image[3,:,:], out_image[2,:,:]))
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.
- Jensen, John R. Remote sensing of the environment: An earth resource perspective 2/e. Pearson Education India, 2009.