Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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!

share|improve this question
add comment

2 Answers

up vote 2 down vote accepted

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

share|improve this answer
    
I knew there must be a better way! Now I can freely use the i.sharpen module. Thanks for pointing out the r.rescale module. Awesome work Nikos! –  Haziq Jul 30 '13 at 7:36
    
Haziq, I am not sure if and how much of the "fine" details are lost when converting 11-bit data sets to 8-bit. QuickBird imagery, for example, is an 11-bit sensor. They are available in both 8-bit and 16-bit formats. It is up to the user to decide what to do. It would certainly be nice for i.pansharpen to handle all kinds of formats. Please have a look at a related "ticket": Ticket #2048: i.pansharpen limited to 8-bit imagery. On the other hand, I simply might not understand stuff and, thus, not able to use i.pansharpen properly... ? –  Nikos Alexandris Jul 30 '13 at 11:10
add comment

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:

Before

After pan-sharpening:

After

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

share|improve this answer
1  
Thanks for your research and efforts, I also had this problem in new GRASS 7 from svn. Now the colours of pansharpened image is OK. =) –  Vladimir Nov 4 '12 at 6:50
    
Oh yeah, I forgot to mention which version of GRASS I'm running - thanks @VladimirNaumov for reminding! I am using GRASS 7svn; should have realized that this issue could be something that doesn't happen across all versions of GRASS (I haven't tried other versions).. –  Haziq Nov 4 '12 at 14:51
    
Great research effort, thank you. –  Nikos Alexandris Jan 4 '13 at 22:59
    
Back to this: you can convert floating point raster maps to integer raster maps by using r.mapcalc along with its integrated int() function. I think r.recode is not meant to be used in the context of your application. –  Nikos Alexandris Jul 28 '13 at 22:13
    
I think the overall "problem" regarding i.pansharpen is that it expects 8-bit raster maps as input, i.e. the input values should range from 0 to 255. –  Nikos Alexandris Jul 28 '13 at 22:46
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.