近年のリモートセンシングデータは、パンクロマティックセンサーによって高解像度のモノクロ画像（＝パンクロマティック画像）を取得できます。一方、解像度は少し劣るものの、マルチスペクトルセンサーからはカラー情報が得られるため、この両者の特性を擬似的に合成させることにより、高解像度カラー画像の作成が可能となります。この方法がパンシャープン処理です。
※合成させるデータ間の空間分解能の差は2倍〜4倍程度が適当とされている。

HSV変換

カラー画像の色情報は、色相（Hue）と彩度（Saturation）と明度（Value）の成分に変換することができると同時に、その逆への交換も可能。
3つのうち、Value成分の画像はパンクロマティック画像と交換して、Hue成分とSaturation成分とパンクロマティック画像からRGBに逆変換することができる。このようにして、擬似的に高解像度のカラー画像を得る手法がHSV変換である。

Google Earth Engineにはすでにee.Image.rgbToHsvが用意されているので、簡単にパンシャープン画像が作成できます。
今回は2020年1月〜2020年12月までのLandsat 8から東京駅周辺の画像を取得してみます。

import ee
import geemap
ee.Initialize()

center_lat = 35.68140174080035; center_lon = 139.76713553334295 # 東京駅

Image = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA') \
 .filterDate('2020-01-01', '2020-12-31').filterBounds(ee.Geometry.Point(center_lon-360, center_lat)).sort('CLOUD_COVER').first())

rgb = image.select('B4', 'B3', 'B2')
pan = image.select('B8')

​sharpened = ee.Image.cat(rgb.rgbToHsv().select('hue', 'saturation'), pan).hsvToRgb()
​
Map = geemap.Map(center=[center_lat, center_lon-360], zoom=14)
left_layer = geemap.ee_tile_layer(rgb, {'max': 0.4}, 'Original')
right_layer = geemap.ee_tile_layer(sharpened, {'max': 0.4}, 'Pansharpened')
Map.split_map(left_layer, right_layer)
Map

合成後の色合いも自然であり、判読用には有効な手法である。

参考になりそうな動画貼っておきます。

Brovey変換

こちらのサイトによると、

この変換手法では、画像中の輝度の差が大きい物体が強調され、建造物のような高輝度部分と影や水域部分とのコントラストが際立って表現されるため、都市域を含む画像に適しています。

といっても、gdal_pansharpen.py（現時点でサポートされているのは “weighted” Brovey algorithmのみらしい）で一発です。

region = ee.Geometry.Rectangle([center_lon-0.1, center_lat-0.1, center_lon+0.1, center_lat+0.1])
yyyymmdd = image.getInfo()['id'][-8:]

# 各バンドをGeoTIFF形式で保存
geemap.ee_export_image(image.select('B8'), f'B8_{yyyymmdd}.tif', crs='epsg:32654', scale=15, region=region)
geemap.ee_export_image(image.select('B2'), f'B2_{yyyymmdd}.tif', crs='epsg:32654', scale=30, region=region)
geemap.ee_export_image(image.select('B3'), f'B3_{yyyymmdd}.tif', crs='epsg:32654', scale=30, region=region)
geemap.ee_export_image(image.select('B4'), f'B4_{yyyymmdd}.tif', crs='epsg:32654', scale=30, region=region)
geemap.ee_export_image(image.select('B5'), f'B5_{yyyymmdd}.tif', crs='epsg:32654', scale=30, region=region)

# パンシャープン処理
import subprocess
cmd = f'gdal_pansharpen.py B8_{yyyymmdd}.tif B4_{yyyymmdd}.tif B2_{yyyymmdd}.tif B3_{yyyymmdd}.tif pan-sharpen.tif'
subprocess.run([cmd], shell=True)

出力したpan-sharpen.tifを描画してみましょう。

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

pansharpen = xr.open_rasterio('pan-sharpen.tif')

red_pan = pansharpen.values[0]
blue_pan = pansharpen.values[1]
green_pan = pansharpen.values[2]

blue = xr.open_rasterio(f'B2_{yyyymmdd}.tif').values[0]
green = xr.open_rasterio(f'B3_{yyyymmdd}.tif').values[0]
red = xr.open_rasterio(f'B4_{yyyymmdd}.tif').values[0]
ni = xr.open_rasterio(f'B5_{yyyymmdd}.tif').values[0]

def convert(band, gamma):
   band = np.clip(band, 0, 1)
   band = np.power(band, 1/gamma)
   return band

red_pan = convert(red_pan, 2.2)
blue_pan = convert(blue_pan, 2.2)
green_pan = convert(green_pan, 2.2)

red = convert(red, 2.2)
blue = convert(blue, 2.2)
green = convert(green, 2.2)

x_min = 250; x_max = 400; y_min = 250; y_max = 400

plt.figure(figsize=(20,20))

plt.subplot(1,2,1)
plt.imshow(np.dstack([red[y_min:y_max, x_min:x_max], green[y_min:y_max, x_min:x_max], blue[y_min:y_max, x_min:x_max]]), vmax=0.4)
plt.title("RGB")

plt.subplot(1,2,2)
plt.imshow(np.dstack([red_pan[y_min*2:y_max*2, x_min*2:x_max*2], green_pan[y_min*2:y_max*2, x_min*2:x_max*2], blue_pan[y_min*2:y_max*2, x_min*2:x_max*2]]), vmax=0.4)
plt.title("PAN RGB")

左側がオリジナルで、右側はパンシャープン画像

↓こちらのnotebookが大変参考になりました！

複数バンドのスペクトル情報を活用したHSV変換とBrovey変換のほか、複数バンドの色の分布（空間情報）を活用した主成分分析法やWavelet変換もあるようです。