This notebook is part of the collection accompanying the talk “Analyzing Satellite Images With Python Scientific Stack” by Milos Miljkovic given at PyData NYC 2014. The content is BSD licensed.
Panchromatic sharpening is a radiometric transformation which uses a higher resolution panchromatic band to fuse it with lower resolution red-green-blue bands. In other words, spatial information in the high-resolution grayscale Band 8 and color information in the RGB color Bands 4, 3, and 2 is merged to create a high-resolution color image.
There are a few algorithms for sharpening: Brovey, IHS, PCA, Gram-Schmidt, simple mean, and weighted algorithms used on per satellite data sets.
In [1]:
# To install watermark extension execute:
# %install_ext https://raw.githubusercontent.com/HyperionAnalytics/watermark/master/watermark.py
%load_ext watermark
In [2]:
%watermark -v -m -p numpy,scikit-image,matplotlib
In [3]:
cd /Users/kronos/gis/l8/guinea_bissau/
In [4]:
import numpy as np
from skimage import io, transform, exposure
from matplotlib import pyplot as plt
In [5]:
%matplotlib inline
In [6]:
def read_band(n):
"""
Load Landsat 8 band
Input:
n - integer in the range 1-11
Output:
img - 2D array of uint16 type
"""
if n in range(1, 12):
tif_list = !ls *.TIF
band_name = 'B' + str(n) + '.TIF'
img_idx = [idx for idx, band_string in enumerate(tif_list) if band_name in band_string]
img = io.imread(tif_list[img_idx[0]])
return img
else:
print('Band number has to be in the range 1-11!')
In [7]:
b4 = read_band(4)
b3 = read_band(3)
b2 = read_band(2)
b8 = read_band(8)
In [8]:
img432_roi = np.dstack((b4, b3, b2))[4400:5200, 3500:4100, :]
In [9]:
b4 = b4/b4.max()
b3 = b3/b3.max()
b2 = b2/b2.max()
b8 = b8/b8.max()
In [10]:
b4 = b4[4400:5200, 3500:4100]
b3 = b3[4400:5200, 3500:4100]
b2 = b2[4400:5200, 3500:4100]
b8 = b8[8800:10400, 7000:8200]
In [11]:
img432 = np.dstack((b4, b3, b2))
img432_2x = transform.rescale(img432, 2)
In [12]:
m = np.sum(img432_2x, axis=2)
ps4 = b8*img432_2x[:, :, 0]/m
ps3 = b8*img432_2x[:, :, 1]/m
ps2 = b8*img432_2x[:, :, 2]/m
img432_ps = np.dstack((ps4, ps3, ps2))
In [13]:
fig = plt.figure(figsize=(12, 9))
fig.set_facecolor('white')
ax1 = fig.add_subplot(121)
ax1.imshow(img432_ps)
plt.title('Pansharpened image')
plt.axis('off')
ax2 = fig.add_subplot(122)
ax2.imshow(img432_roi/65535)
plt.title('Raw image')
plt.axis('off')
plt.show()
In [14]:
fig = plt.figure(figsize=(10, 7))
fig.set_facecolor('white')
for color, channel in zip('rgb', np.rollaxis(img432_ps, axis=-1)):
counts, centers = exposure.histogram(channel)
plt.plot(centers[1::], counts[1::], color=color)
plt.show()
In [15]:
fig = plt.figure(figsize=(10, 7))
fig.set_facecolor('white')
for color, channel in zip('rgb', np.rollaxis(img432_roi, axis=-1)):
counts, centers = exposure.histogram(channel)
plt.plot(centers[1::], counts[1::], color=color)
plt.show()
In [16]:
img1 = np.empty(img432_ps.shape, dtype='float64')
lims = [(0.088,0.17), (0.108, 0.19), (0.130,0.20)]
for lim, channel in zip(lims, range(3)):
img1[:, :, channel] = exposure.rescale_intensity(img432_ps[:, :, channel], lim)
In [17]:
img2 = np.empty(img432_roi.shape, dtype='float64')
lims = [(7100,14500), (8200, 14000), (9200,13500)]
for lim, channel in zip(lims, range(3)):
img2[:, :, channel] = exposure.rescale_intensity(img432_roi[:, :, channel], lim)
In [18]:
fig = plt.figure(figsize=(12, 9))
fig.set_facecolor('white')
ax1 = fig.add_subplot(121)
ax1.imshow(img1)
plt.title('Pansharpened image')
plt.axis('off')
ax2 = fig.add_subplot(122)
ax2.imshow(img2/65535)
plt.title('Raw image')
plt.axis('off')
plt.show()
In [19]:
fig = plt.figure(figsize=(9, 12))
fig.set_facecolor('white')
ax1 = fig.add_subplot(121)
ax1.imshow(img1[750:1050, 880:1040, :])
plt.title('Pansharpened image - zoomed')
plt.axis('off')
ax2 = fig.add_subplot(122)
ax2.imshow(img2[375:525, 440:520, :]/65535)
plt.title('Raw image - zoomed')
plt.axis('off')
plt.show()
This notebook is part of the collection accompanying the talk “Analyzing Satellite Images With Python Scientific Stack” by Milos Miljkovic given at PyData NYC 2014. The content is BSD licensed.
Panchromatic sharpening is a radiometric transformation which uses a higher resolution panchromatic band to fuse it with lower resolution red-green-blue bands. In other words, spatial information in the high-resolution grayscale Band 8 and color information in the RGB color Bands 4, 3, and 2 is merged to create a high-resolution color image.
There are a few algorithms for sharpening: Brovey, IHS, PCA, Gram-Schmidt, simple mean, and weighted algorithms used on per satellite data sets.
In [1]:
# To install watermark extension execute:
# %install_ext https://raw.githubusercontent.com/HyperionAnalytics/watermark/master/watermark.py
%load_ext watermark
In [2]:
%watermark -v -m -p numpy,scikit-image,matplotlib
In [3]:
cd /Users/kronos/gis/l8/guinea_bissau/
In [4]:
import numpy as np
from skimage import io, transform, exposure
from matplotlib import pyplot as plt
In [5]:
%matplotlib inline
In [6]:
def read_band(n):
"""
Load Landsat 8 band
Input:
n - integer in the range 1-11
Output:
img - 2D array of uint16 type
"""
if n in range(1, 12):
tif_list = !ls *.TIF
band_name = 'B' + str(n) + '.TIF'
img_idx = [idx for idx, band_string in enumerate(tif_list) if band_name in band_string]
img = io.imread(tif_list[img_idx[0]])
return img
else:
print('Band number has to be in the range 1-11!')
In [7]:
b4 = read_band(4)
b3 = read_band(3)
b2 = read_band(2)
b8 = read_band(8)
In [8]:
img432_roi = np.dstack((b4, b3, b2))[4400:5200, 3500:4100, :]
In [9]:
b4 = b4/b4.max()
b3 = b3/b3.max()
b2 = b2/b2.max()
b8 = b8/b8.max()
In [10]:
b4 = b4[4400:5200, 3500:4100]
b3 = b3[4400:5200, 3500:4100]
b2 = b2[4400:5200, 3500:4100]
b8 = b8[8800:10400, 7000:8200]
In [11]:
img432 = np.dstack((b4, b3, b2))
img432_2x = transform.rescale(img432, 2)
In [12]:
m = np.sum(img432_2x, axis=2)
ps4 = b8*img432_2x[:, :, 0]/m
ps3 = b8*img432_2x[:, :, 1]/m
ps2 = b8*img432_2x[:, :, 2]/m
img432_ps = np.dstack((ps4, ps3, ps2))
In [13]:
fig = plt.figure(figsize=(12, 9))
fig.set_facecolor('white')
ax1 = fig.add_subplot(121)
ax1.imshow(img432_ps)
plt.title('Pansharpened image')
plt.axis('off')
ax2 = fig.add_subplot(122)
ax2.imshow(img432_roi/65535)
plt.title('Raw image')
plt.axis('off')
plt.show()
In [14]:
fig = plt.figure(figsize=(10, 7))
fig.set_facecolor('white')
for color, channel in zip('rgb', np.rollaxis(img432_ps, axis=-1)):
counts, centers = exposure.histogram(channel)
plt.plot(centers[1::], counts[1::], color=color)
plt.show()
In [15]:
fig = plt.figure(figsize=(10, 7))
fig.set_facecolor('white')
for color, channel in zip('rgb', np.rollaxis(img432_roi, axis=-1)):
counts, centers = exposure.histogram(channel)
plt.plot(centers[1::], counts[1::], color=color)
plt.show()
In [16]:
img1 = np.empty(img432_ps.shape, dtype='float64')
lims = [(0.088,0.17), (0.108, 0.19), (0.130,0.20)]
for lim, channel in zip(lims, range(3)):
img1[:, :, channel] = exposure.rescale_intensity(img432_ps[:, :, channel], lim)
In [17]:
img2 = np.empty(img432_roi.shape, dtype='float64')
lims = [(7100,14500), (8200, 14000), (9200,13500)]
for lim, channel in zip(lims, range(3)):
img2[:, :, channel] = exposure.rescale_intensity(img432_roi[:, :, channel], lim)
In [18]:
fig = plt.figure(figsize=(12, 9))
fig.set_facecolor('white')
ax1 = fig.add_subplot(121)
ax1.imshow(img1)
plt.title('Pansharpened image')
plt.axis('off')
ax2 = fig.add_subplot(122)
ax2.imshow(img2/65535)
plt.title('Raw image')
plt.axis('off')
plt.show()
In [19]:
fig = plt.figure(figsize=(9, 12))
fig.set_facecolor('white')
ax1 = fig.add_subplot(121)
ax1.imshow(img1[750:1050, 880:1040, :])
plt.title('Pansharpened image - zoomed')
plt.axis('off')
ax2 = fig.add_subplot(122)
ax2.imshow(img2[375:525, 440:520, :]/65535)
plt.title('Raw image - zoomed')
plt.axis('off')
plt.show()