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

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..