Pansharpening

Pansharpening#

Landsat includes a panchromatic band which has greater spatial resolution (15 m versus the standard 30 m of visible bands). It can be used to artificially increase the resolution of the visible bands (red, green, and blue) in a process called pansharpening.

To illustrate this concept, let’s download a scene from just North of the city of Liverpool, UK:

import xlandsat as xls
import matplotlib.pyplot as plt

path = xls.datasets.fetch_liverpool()
path_pan = xls.datasets.fetch_liverpool_panchromatic()

Load the scene with xlandsat.load_scene:

scene = xls.load_scene(path)
scene
<xarray.Dataset> Size: 2MB
Dimensions:          (easting: 433, northing: 267)
Coordinates:
  * easting          (easting) float64 3kB 4.87e+05 4.87e+05 ... 5e+05 5e+05
  * northing         (northing) float64 2kB 5.922e+06 5.922e+06 ... 5.93e+06
Data variables:
    coastal_aerosol  (northing, easting) float16 231kB 0.06238 ... 0.0769
    blue             (northing, easting) float16 231kB 0.0708 ... 0.08301
    green            (northing, easting) float16 231kB 0.08618 ... 0.09753
    red              (northing, easting) float16 231kB 0.06824 0.06824 ... 0.116
    nir              (northing, easting) float16 231kB 0.04553 0.04553 ... 0.173
    swir1            (northing, easting) float16 231kB 0.04626 ... 0.2213
    swir2            (northing, easting) float16 231kB 0.047 0.04663 ... 0.1686
    thermal          (northing, easting) float16 231kB 287.0 287.0 ... 290.5
Attributes: (12/19)
    Conventions:                CF-1.8
    title:                      Landsat 8 scene from 2020-09-27 (path/row=204...
    digital_object_identifier:  https://doi.org/10.5066/P9OGBGM6
    origin:                     Image courtesy of the U.S. Geological Survey
    landsat_product_id:         LC08_L2SP_204023_20200927_20201006_02_T1
    processing_level:           L2SP
    ...                         ...
    ellipsoid:                  WGS84
    date_acquired:              2020-09-27
    scene_center_time:          11:10:50.3140030Z
    wrs_path:                   204
    wrs_row:                    23
    mtl_file:                   GROUP = LANDSAT_METADATA_FILE\n  GROUP = PROD...

And load the panchromatic band with xlandsat.load_panchromatic:

panchromatic = xls.load_panchromatic(path_pan)
panchromatic
<xarray.DataArray 'pan' (northing: 534, easting: 867)> Size: 2MB
array([[0.04228   , 0.04231999, 0.04229999, ..., 0.0603    , 0.0497    ,
        0.04848   ],
       [0.04224   , 0.04207999, 0.04212   , ..., 0.05284   , 0.04792   ,
        0.03968   ],
       [0.04251999, 0.04236   , 0.04207999, ..., 0.05194   , 0.04589999,
        0.03664   ],
       ...,
       [0.04477999, 0.04421999, 0.04423999, ..., 0.05402   , 0.05444   ,
        0.05415999],
       [0.04465999, 0.04457999, 0.04387999, ..., 0.05435999, 0.05398   ,
        0.05471999],
       [0.04421999, 0.04408   , 0.04382   , ..., 0.05731999, 0.05505999,
        0.05444   ]], shape=(534, 867), dtype=float32)
Coordinates:
  * easting   (easting) float64 7kB 4.87e+05 4.87e+05 4.87e+05 ... 5e+05 5e+05
  * northing  (northing) float64 4kB 5.922e+06 5.922e+06 ... 5.93e+06 5.93e+06
Attributes: (12/25)
    long_name:                  panchromatic
    units:                      reflectance
    number:                     8
    filename:                   LC08_L1TP_204023_20200927_20201006_02_T1_B8.TIF
    scaling_mult:               2e-05
    scaling_add:                -0.1
    ...                         ...
    ellipsoid:                  WGS84
    date_acquired:              2020-09-27
    scene_center_time:          11:10:50.3140030Z
    wrs_path:                   204
    wrs_row:                    23
    mtl_file:                   GROUP = LANDSAT_METADATA_FILE\n  GROUP = PROD...

Now we can plot an RGB composite and the panchromatic band for comparison:

rgb = xls.composite(scene, rescale_to=(0, 0.25))

plt.figure(figsize=(16, 10))
ax = plt.subplot(2, 1, 1)
rgb.plot.imshow(ax=ax)
ax.set_aspect("equal")
ax.set_title("RGB")
ax = plt.subplot(2, 1, 2)
panchromatic.plot.pcolormesh(
    ax=ax, cmap="gray", vmin=0.02, vmax=0.1, add_colorbar=False,
)
ax.set_aspect("equal")
ax.set_title("Panchromatic")
plt.tight_layout()
_images/pansharpen_3_0.png

The pansharpening is implemented in xlandsat.pansharpen:

scene_sharp = xls.pansharpen(scene, panchromatic)
scene_sharp
<xarray.Dataset> Size: 6MB
Dimensions:   (northing: 534, easting: 867)
Coordinates:
  * northing  (northing) float64 4kB 5.922e+06 5.922e+06 ... 5.93e+06 5.93e+06
  * easting   (easting) float64 7kB 4.87e+05 4.87e+05 4.87e+05 ... 5e+05 5e+05
Data variables:
    red       (northing, easting) float32 2MB 0.03763 0.03767 ... 0.06035
    green     (northing, easting) float32 2MB 0.04753 0.04758 ... 0.05076
    blue      (northing, easting) float32 2MB 0.03905 0.03909 ... 0.04369 0.0432
Attributes: (12/21)
    Conventions:                  CF-1.8
    title:                        Pansharpend Landsat 8 scene from 2020-09-27...
    digital_object_identifier:    https://doi.org/10.5066/P9OGBGM6
    origin:                       Image courtesy of the U.S. Geological Survey
    landsat_product_id:           LC08_L2SP_204023_20200927_20201006_02_T1
    processing_level:             L2SP
    ...                           ...
    scene_center_time:            11:10:50.3140030Z
    wrs_path:                     204
    wrs_row:                      23
    pansharpening_method:         Weighted Brovey Transform
    pansharpening_rgb_weights:    (1, 1, 0.2)
    pansharpening_band_filename:  LC08_L1TP_204023_20200927_20201006_02_T1_B8...

Finally, let’s compare the sharpened and original RGB composites:

rgb_sharp = xls.composite(scene_sharp, rescale_to=(0, 0.15))

plt.figure(figsize=(16, 10))
ax = plt.subplot(2, 1, 1)
rgb.plot.imshow(ax=ax)
ax.set_aspect("equal")
ax.set_title("Original")
ax = plt.subplot(2, 1, 2)
rgb_sharp.plot.imshow(ax=ax)
ax.set_aspect("equal")
ax.set_title("Pansharpened")
plt.tight_layout()
_images/pansharpen_5_0.png