How to color-balance false-composites of pan-sharpened Landsat images in GRASS?

Question

Whenever I try to pan-sharpen composites of some Landsat images in GRASS using i.pansharpen, i.fusion.brovey or the IHS sharpening method, the output will have some or all of the following characteristics:

• the composite color is in a different hue compared to the un-sharpened composite
• the brightness level is messed up
• the entire composite went all-white/all-black (when using images pre-processed to top-of-atmosphere reflectance or surface reflectance corrections in i.landsat.toar)

I've also tried all of the following; but the colors/brightness remained the same or turned even worse:

• Applied i.landsat.rgb, before-and-after the pan-sharpening process
• Played with the -f or -p flag in i.landsat.rgb
• Tried r.colors to edit the color table to grey/grey255/grey.eq
• Tried i.pansharpen using all Brovey/IHS/PCA methods
• Played with the -l flag in i.pansharpen to rebalance the blue-channel

The GRASS GIS manual did explained on how to perform pan-sharpening and color-balancing, but I can't figure out how to combine both processes in a concurrent workflow. I suspected that this is due to my poor understanding of color-tables, color-histogram, etc. in GRASS..

So, can someone explain to me - how do you tackle color-balancing problems when dealing with Landsat images after image-processing in GRASS? Can you share with me your favorite workflow/methods?

Many thanks for any feedback!

Answer 1

Overview

One working approach inside GRASS-GIS version 7 to get an acceptable color-balanced composite image after Pan-sharpening is

1. check if input data are 8-bit ranging inside [0, 255]
2. if the data are inside [0, 255] proceed then to pan-sharpening (i.pansharpen)
3. if the data are not inside [0, 255], rescale them to this range (r.rescale)
4. pan-sharpen with any of the featured methods (Brovey, IHS, PCA)
5. color-balance automatically by using the i.landsat.rgb module or manually adjusting the color tables of the bands of interest

Details and example instructions

Pan-Sharpening / Fusion

GRASS 7 holds a dedicated pan-sharpening module, i.pansharpen which features three techniques for sharpening, namely the Brovey transformation, the classical IHS method and one that is based on PCA.

i.pansharpen works fine with 8-bit raster maps as an input. If the data to be processed are out of this range, that is out of [0, 255], they can be rescaled to fit into this range by using GRASS' r.rescale module.

Given a set of 11-bit spectral bands (for example Blue, Green, Red, NIR and Pan) ranging between [0, 2047], querying the Blue band for example would return

r.info Blue_DNs -r
min=0
max=2047

Rescaling the Blue band to range between [0, 255]

r.rescale in=Blue_DNs out=Blue_DNs_255 from=0,2047 to=0,255

The same step applies to both the rest of the multi-spectral bands and the Panchromatic band of interest.

As usual when working with GRASS, it is required to set the region of interest, i.e. g.regionrast=Blue_DNs_255 to match the extent of the band(s) or else. The resolution itself is taken care in this particular case by the module and the resulting pan-sharpened raster maps will be of the same high(er) resolution as the Panchromatic band.

An example command for an IHS-based Pan-Sharpening action might look like

i.pansharpen pan=Pan_DNs_255 ms1=Blue_DNs_255 ms2=Green_DNs_255 ms3=Red_DNs_255 output=sharptest255 sharpen=ihs

Color Balancing

After the process completion, the module outputs

...
The following pan-sharpened output maps have been generated:
sharptest255_red
sharptest255_green
sharptest255_blue

To visualize output, run: g.region -p rast=sharptest255.red
d.rgb r=sharptest255_red g=sharptest255_green b=sharptest255_blue

Normally it should be enough to re-balance the colors after the pan-sharpening by using for example the i.landsat.rgb module or manual adjustment of each of the three bands that would compose an RGB image.

Screenshots

...to be added

Answer 2

I've searched high-and-low, and I think I've discovered the root of my problems. I believe I got the solution for them now - but it's a little bit messy. I'm sure there are better ways to solve them. Do share if you know an easier way!

ROOTS OF PROBLEMS:

1. The output of i.landsat.toar is in floating point. I've realized that when I use floating point rasters in any pan-sharpening method, the colors will mess-up. Those algorithm somehow preferred rasters in the original integer form.
2. Pan-sharpening modules such as i.pansharpen and i.fusion.brovey modules will mess-up the colors. I haven't quite grasp the algorithms that they used in those modules - but somehow the color-tables will be affected, and ruining the resulting pan-sharpened images.

SOLUTION:

1. Convert output from i.landsat.toar from float to int, using r.recode.
2. Use the rasters as inputs in IHS pan-sharpening method (i.rgb.his and i.his.rgb). I keep away from using i.pansharpen and i.fusion.brovey.

THE WORKFLOW:

1. Use r.info with the -r flag to get the DN min and max values of each raster bands which have been processed with i.landsat.toar. For example:

   > r.info -r BAND1
   min=0.01
   max=0.370064120902708
   

   As we can see, the values is between 0-1, which are pretty different than the original ones (which are between 0-255). That explains why the output from pan-sharpening turned out blank, because the used range-of-value is very low (below 1).

2. Convert that raster band using r.recode. Use the min and max values obtained from step 1 to convert into a new range of 0-255. An example code snippet:

   r.recode input=BAND1 output=NEWBAND1 rules=- << EOF
   0.01:0.370064120902708:0:255
   EOF
   

   We can check the new converted values with r.info:

   > r.info -r NEWBAND1
   min=1
   max=254
   

   Values are in 0-255: now it's usable for pan-sharpening process.

3. Apply gray-scale color table to the converted band with r.colors.

   r.colors NEWBAND1 color=grey
   

   So far, I get the best results using the grey color table - the pan-sharpened composites matched closely with the original composites. The other alternatives are to equalize grey color table with color=grey.eq or using the -e flag with color=grey. Or we can use the i.landsat.rgb module instead of r.colors..

4. Repeat step 1-3 with other raster bands that we intend to be use as composites, including the pan raster (band 8). Use of scripts would be much appreciated here.

5. Then use the processed rasters as inputs in IHS pan-sharpening method. For example, when making the composite of band 7,4,2:

   i.rgb.his r=NEWBAND7 g=NEWBAND4 b=NEWBAND2 hue=HUE int=INT sat=SAT
   

   This will output 3 layers: a hue layer HUE, an intensity layer INT, and also a saturation layer SAT. We will then replace the intensity layer INT with the pan raster band NEWBAND8 in i.his.rgb:

   i.his.rgb hue=HUE sat=SAT int=NEWBAND8 r=COMP742_red g=COMP742_green b=COMP742_blue
   

   Resulting red channels of COMP742_red, COMP742_green, COMP742_blue can then be combined using d.rgb or r.composite..

SAMPLE BEFORE & AFTER:

Before pan-sharpening:

After pan-sharpening:

Maybe it's hard to tell the sharpening differences, when viewing from such small images. But what's important is the color of the pan-sharpened image matched the composite from the original. Mission accomplished!

OTHER NOTES:

• Don't r.recode the thermal bands (band 6). i.landsat.toar output these thermal bands in Kelvin temperature values (nothing to do with DN values). Keep the r.recode routine on the normal multi-spectral and pan bands (bands 1-5,7,8).
• If we never even use i.landsat.toar, but the resulting composites look really wrong, it's usually because of mismatching of color-tables before and after pan-sharpening process. I applied r.colors RASTER color=grey to the original raster bands before pan-sharpening, and to the resulting channels after pan-sharpening to ensure close-matching of colors.
• A usual case of the wrong-color composite problem: the original raster bands are in color=grey255, the output of i.landsat.rgb is in color=grey.eq. No wonder they both look different!
• Processing landsat images for use could really be a time-taxing activity. Better have something to do when waiting for it all to process, or at least have some ludicrous amount of coffee and good music while you're at it ;)

Hope this would benefit someone : took me days to find what's wrong..