pansharpening

[aloboa]

Could pansharpening (eg. the one in gdal_pansharpen.py) be implemented?

[rhijmans]

I want to give it a try. Could you provide a small example input data set with expected output?

[kadyb]

Here you can download sample Landsat data without registration. Band 8 is a panchromatic band in Landsat 7 and 8. Below I wrote a simple resampling function using the panchromatic band.

# Landsat 7 sample data
# https://landsat.usgs.gov/sites/default/files/C2_Sample_Data/LE07_L1TP_042027_20050927_20200914_02_T1.zip

library("terra")

path = "LE07_L1TP_042027_20050927_20200914_02_T1/"
imgs = list.files(path, ".+B[1234578]+\\.TIF$", full.names = TRUE)

blue = rast(imgs[1])
green = rast(imgs[2])
red = rast(imgs[3])
pan = rast(imgs[7])

pansharpen = function(B, G, R, PAN) {

  B_hires = resample(B, PAN, method = "near")
  G_hires = resample(G, PAN, method = "near")
  R_hires = resample(R, PAN, method = "near")

  calc = function(a, b) (a + b) / 2
  B_out = lapp(c(B_hires, PAN), fun = calc)
  G_out = lapp(c(G_hires, PAN), fun = calc)
  R_out = lapp(c(R_hires, PAN), fun = calc)

  output = c(B_out, G_out, R_out)
  names(output) = c("B", "G", "R")
  return(output)

}

RGB_30m = c(blue, green, red)
RGB_15m = pansharpen(blue, green, red, pan)

res(RGB_30m)
#> [1] 30 30
res(RGB_15m)
#> [1] 15 15

There are also available some functions in RStoolbox and satelliteTools packages.

[aloboa]

Using a Landsat subscene is a good start, I will prepare and provide a hyperspectral example also in a short time.
@kadyb Pansharpening is not just an interpolation as done by resample(), it is an interpolation guided by another image, most often a panchromatic image.

The RStools function implements the most traditional approaches (pca, Brovey, IHS), but still uses raster instead of terra. Note IHS is valid for panchromatic + RGB imagery only.
According to
https://gdal.org/drivers/raster/vrt.html#gdal-vrttut-pansharpen
gdal_pansharpen.py applies the Brovey algorithm, but have not tested the actual differences with RStolls::panSharpen
The satellite tools function is more sophisticated (but have not tested if any better than the traditional methods), implementing OTB methods (RCS (==Brovey), LMVM, Bayesian). It looks like this package is no longer maintained.

I find particularly enlightening the following:
Mhangara, Paidamwoyo, Willard Mapurisa, and Naledzani Mudau. “Comparison of Image Fusion Techniques Using Satellite Pour l’Observation de La Terre (SPOT) 6 Satellite Imagery.” Applied Sciences 10, no. 5 (March 10, 2020): 1881. https://doi.org/10.3390/app10051881.
Vivone, Gemine, Luciano Alparone, Jocelyn Chanussot, Mauro Dalla Mura, Andrea Garzelli, Giorgio A. Licciardi, Rocco Restaino, and Lucien Wald. “A Critical Comparison Among Pansharpening Algorithms.” IEEE Transactions on Geoscience and Remote Sensing 53, no. 5 (May 2015): 2565–86. https://doi.org/10.1109/TGRS.2014.2361734.
Loncan, Laetitia, Luis B. Almeida, José M. Bioucas-Dias, Xavier Briottet, Jocelyn Chanussot, Nicolas Dobigeon, Sophie Fabre, et al. “Hyperspectral Pansharpening: A Review.” ArXiv:1504.04531 [Physics, Stat], April 17, 2015. http://arxiv.org/abs/1504.04531.

[kadyb]

@kadyb Pansharpening is not just an interpolation as done by resample(), it is an interpolation guided by another image, most often a panchromatic image.

Note that in my example, not only resampling is performed, but averaging the values with the panchromatic band is also included. You can change it to a Brovey formula using this code (but it would probably be better to include weights to limit the influence of the blue band):

calc = function(a, b, c, d) a / (a + b + c) * d
B_out = lapp(c(B_hires, G_hires, R_hires, PAN), fun = calc)
G_out = lapp(c(G_hires, B_hires, R_hires, PAN), fun = calc)
R_out = lapp(c(R_hires, B_hires, G_hires, PAN), fun = calc)

but still uses raster instead of terra

At this stage, I wouldn't see it as a disadvantage. {raster} is faster than {terra} in some cases (e.g. raster arithmetic) and benchmarks can help here. Another question is the {terra} package a good place for such advanced tools for the pansharpening or should they be in a separate package?

[aloboa]

Do you mean terra is not the place for advanced methods? I disagree. Pansharpening is a very specific operation, if we
split tasks into different packages at this level we would end up with a very inconvenient (from a users's point of view) fragmentation.
Anyway, this is obviously up to the developer.

[rhijmans]

I have now implemented two algorithms in luna::pansharpen. Depending on how develops that may be a good place for it, but otherwise it can go back to terra. I think I can speed things up a bit more, but here is a start.

[aloboa]

Apologies for my late answer.
I have tested and compared to RStoolbox::panSharpen() and gdal_pansharpen.py
https://www.dropbox.com/s/6lti3n5vhf1fexz/pansharpout.pdf?dl=0
I think that the bias vs. the original is caused by the calculated low-resolution panchromatic image being biased vs. the actual input high-res panchromatic image. Could you let the user input the weights? Also, could you implement the PCA method?