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>
Dimensions:          (easting: 433, northing: 267)
Coordinates:
  * easting          (easting) float64 4.87e+05 4.87e+05 ... 5e+05 5e+05
  * northing         (northing) float64 5.922e+06 5.922e+06 ... 5.93e+06
Data variables:
    coastal_aerosol  (northing, easting) float16 0.06238 0.06238 ... 0.0769
    blue             (northing, easting) float16 0.0708 0.07117 ... 0.08301
    green            (northing, easting) float16 0.08618 0.08667 ... 0.09753
    red              (northing, easting) float16 0.06824 0.06824 ... 0.116
    nir              (northing, easting) float16 0.04553 0.04553 ... 0.173
    swir1            (northing, easting) float16 0.04626 0.04651 ... 0.2213
    swir2            (northing, easting) float16 0.047 0.04663 ... 0.1744 0.1686
    thermal          (northing, easting) float16 287.0 287.0 ... 290.5 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)>
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   ]], dtype=float32)
Coordinates:
  * easting   (easting) float64 4.87e+05 4.87e+05 4.87e+05 ... 5e+05 5e+05 5e+05
  * northing  (northing) float64 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>
Dimensions:   (northing: 534, easting: 867)
Coordinates:
  * northing  (northing) float64 5.922e+06 5.922e+06 ... 5.93e+06 5.93e+06
  * easting   (easting) float64 4.87e+05 4.87e+05 4.87e+05 ... 5e+05 5e+05 5e+05
Data variables:
    red       (northing, easting) float32 0.03763 0.03767 ... 0.06103 0.06035
    green     (northing, easting) float32 0.04753 0.04758 ... 0.05133 0.05076
    blue      (northing, easting) float32 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