Preface

In the past 50 years, Landsat series satellites have provided us with very long time series of earth surface observation information. At this stage, Landsat satellites are still in service, providing very valuable data for global governance and scientific research.

Now is the era of big data. As the field of earth science, remote sensing data is absolutely valuable first-hand measured data, and the amount of data is very large. At this stage, the production speed of most remote sensing data may be far faster than that of our use. Therefore, most remote sensing data are still in a sleeping state, waiting for the continuous innovation of technology to wake up the sleeping data and make it glow.

Today, let's explore how to use Landsat data to create a remote sensing image. The picture shows Hangzhou Bay.

The southeast wins, the three Wu capitals, and Qiantang has been prosperous since ancient times.

1.Preparation

First prepare the downloaded Landsat8 data set. Landsat8 includes 11 bands. In the data folder, each band has a tif file. This tif file contains coordinate information, band value and other information, which can be read through the rasterio module of Python (if necessary, you can write about using this library to read the tif file information in detail).

Import the Python library to be used next.

import numpy as np
import os
import matplotlib.pyplot as plt
import tifffile as tiff
from skimage import exposure

When using remote sensing image to make a base map, it is usually to select three bands of data as the three channels of image RGB, where RGB represents three channels of red, green and blue, and different band combinations will get different effects.

We use Landsat data in Google Earth, which has a band combination of 432. Namely, the B4, B3 and B2 bands of Landsat are used as the three channels of the image.

The map of Hangzhou Bay above uses 753. If the band combination is 432, the effect is as shown in the following figure. The true color picture synthesized at 432 band is close to the real color of the ground object, the image is flat and the tone is dark.

2.Start

First, read out the data of each band of Landsat.

l8path = r'pathOfLandsat8\LC08_L1TP_118039_20160126_20170330_01_T1'
fileName = l8path.split("\")[-1]
prefix = os.path.join(l8path,fileName)
bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11']
n_bands = len(bands)
temp = tiff.imread(prefix + "_B1.TIF")
arr = np.zeros((n_bands, np.shape(temp)[0], np.shape(temp)[1]), dtype=np.float)
for i in range(n_bands):
    file2read = prefix + "_" + bands[i] + ".TIF"
    arr[i, :, :] = tiff.imread(file2read)

Select three bands

rgb=(6, 4, 2) #The lower bound of Python is 0, which is actually 753

rgb_bands = arr[rgb, :, :]
rgb_bands = _stretch_im(rgb_bands, 2)
rgb_bands = bytescale(rgb_bands).transpose([1, 2, 0])
plt.close('all')
fig, ax = plt.subplots(figsize=figSize)
ax.imshow(rgb_bands, extent=extent)
ax.set_title(title)
ax.set(xticks=[], yticks=[])

Two functions are used, one is to improve the contrast of the image, and the other is to normalize the band data to the [0255] range commonly used by pixel points.

def _stretch_im(arr, str_clip):
  s_min = str_clip
  s_max = 100 - str_clip
  arr_rescaled = np.zeros_like(arr)
  pLow, pHigh = np.percentile(arr[arr > 0], (s_min, s_max))
  arr_rescaled = exposure.rescale_intensity(arr, in_range=(pLow, pHigh))
  return arr_rescaled.copy()

def bytescale(data, high=255, low=0, cmin=None, cmax=None):
  if data.dtype == "uint8":
    return data
  if high > 255:
    raise ValueError("high should be less than or equal to 255.")
  if low < 0:
    raise ValueError("low should be greater than or equal to 0.")
  if high < low:
    raise ValueError("high should be greater than or equal to low.")
  if cmin is None or (cmin < data.min()):
    cmin = float(data.min())
  if (cmax is None) or (cmax > data.max()):
    cmax = float(data.max())
  # Calculate range of values
  crange = cmax - cmin
  if crange < 0:
    raise ValueError("cmax should be larger than cmin.")
  elif crange == 0:
    raise ValueError(
      "cmax and cmin should not be the same value. Please specify "
      "cmax > cmin"
    )
  scale = float(high - low) / crange
  # If cmax is less than the data max, then this scale parameter will create
  # data > 1.0. clip the data to cmax first.
  data[data > cmax] = cmax
  bytedata = (data - cmin) * scale + low
  return (bytedata.clip(low, high) + 0.5).astype("uint8")

3.Summary

This remote sensing image map can be used as the base map of our map, and can also realize the fusion of two remote sensing image maps, and overlay the data you want to overlay on the image map.