How to use histogram enhancement in yaml files?

[haiyangdaozhang]

Dear:

I checked the satpy documentation and found that it supports histogram enhancement. But I did not find an example or introduction. How can I use histogram enhancement in yaml files or in python codes?

url:
https://satpy.readthedocs.io/en/stable/enhancements.html?highlight=Histogram#histogram

Thank you.

Haidao

[djhoese]

Sorry about that. This should work (based on the linked to documentation for crude stretch):

- name: hist_eq
  method: !!python/name:satpy.enhancements.stretch
  kwargs:
    stretch: histogram

This ends up calling this method in the trollimage package: https://github.com/pytroll/trollimage/blob/5aeef5a2534442434cb00b0a9e091bf72a222130/trollimage/xrimage.py#L1234-L1244

Which takes an additional "approximate" boolean argument which when True will perform a faster approximate histogram equalization. It defaults to False. You can provide this in the YAML at the same level as the stretch argument; so approximate: true should work if you'd like to use it.

[haiyangdaozhang]

I tried to use this enhancement for a compositors, but it seems that it uses histogram enhancement for the entire picture. How do I use it for each channel? Is it can be created in my yaml file？

yam file:

enhancements:
  test:
    name: test
    operations:
    - name: hist_eq
      method: !!python/name:satpy.enhancements.stretch
      kwargs:
        stretch: histogram

[djhoese]

You're correct, the histogram enhancement applies to the entire data array. To apply it per channel you would likely have to write your own compositor python code to apply the enhancement inside the compositor and then make an RGB from that. It isn't incredibly difficult if you are comfortable with python, but it does require more customization than you're looking for. I'm done something kind of similar as far as pre-enhancing a single composite but this was in my own project "polar2grid":

https://github.com/ssec/polar2grid/blob/74c3d56ceb96d2e4890e204bb29fe438a45cf4a7/polar2grid/composites/enhanced.py#L27-L49

I'm not sure how useful that would be for you.

Could you explain what type of product you are making? Wouldn't applying histogram equalization to an RGB make it non-physical? I mean, the colors in the RGB don't actually mean much from image to image. Mostly I'm just curious.

[haiyangdaozhang]

I just want to make the pictures more visible. I use modis channel 1, 2, and 31 , and each channel is histogram equalized to draw smoke.

[djhoese]

In that case, my original suggestions still stand (custom compositor or custom enhancement python code). That, or you define a "crude" or "linear" stretch in your enhancement for your RGB which allows for per-channel min and max limits.

[haiyangdaozhang]

Hello David:

Sorry to disturb you on the weekend.

I tried to use the method you recommended, but I don't know how to use SingleEnhancedBandCompositor.
I guess it should be like the usage of GenericCompositor. Then I tried this, but failed.

comp = SingleEnhancedBandCompositor("test1")
test1 = comp(scn["1"])

When I tried to use SingleBandCompositor, it alse failed

from satpy.composites import SingleBandCompositor
comp = SingleEnhancedBandCompositor("test1")
test1 = comp(scn["1"])

May you show me an example about SingleEnhancedBandCompositor?

Then if I get a SingleEnhancedBandCompositor which is named test1, how can I enhance it ? May be like this ?

test1.gamma(2)

or

test1.stretch_hist_equalize(True)

Thank you for responce.

Haidao

[djhoese]

"It failed" how? Do you have error messages? The output isn't what you expect? How did it fail?

I assumed you wanted to do this all in YAML, sorry. If you're ok doing this all in python then we can probably make something up without having to use the extra classes. It isn't as pretty but it can be done all in one python script so that's nice.

[haiyangdaozhang]

I want to use satpy to create a custom three-channel compositor, and each channel is histogram equalized. You told me that you can use SingleEnhancedBandCompositor.
I guess the script should be like this:

from satpy.writers import to_image
from enhanced import SingleEnhancedBandCompositor
from satpy.composites import GenericCompositor

comp_single = SingleEnhancedBandCompositor("r_band")
r_band = comp_single(scn["1"])
r_band.stretch_hist_equalize(True)

comp_single = SingleEnhancedBandCompositor("g_band")
g_band = comp(scn["2"])
g_band.stretch_hist_equalize(True)

comp_single = SingleEnhancedBandCompositor("b_band")
b_band = comp_single(scn["3"])
b_band.stretch_hist_equalize(True)

comp = GenericCompositor("test1")
test1 = comp((r_band, g_band，b_band))
img = to_image(test1)
img.save("test1.tif")

But when I used the script, I failed at this step.

comp_single = SingleEnhancedBandCompositor("r_band")
r_band = comp_single(scn["1"])

Error messages:

Traceback (most recent call last):
  File "H:\study_program\anaconda3\envs\p2g\lib\site-packages\IPython\core\interactiveshell.py", line 3444, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-85-ff421e32a912>", line 2, in <module>
    r_band = comp_single(scn["1"])
  File "H:\study_program\AAmy_program\python\P2GViewer\program_v4\enhanced.py", line 41, in __call__
    raise ValueError("Can't have more than one band in a single-band composite")
ValueError: Can't have more than one band in a single-band composite

I guess I used SingleEnhancedBandCompositor incorrectly. Then I want to know if the other parts of the script are correct.

Thank you for your help these days.

[djhoese]

So my original suggestion for using SingleEnhancedBandCompositor assumed you wanted to do everything in YAML. But it also requires not only adding things to a composite YAML, but also the enhancement YAML files. So it is pretty involved. Since you're willing to do everything in python this actually can be pretty easy; at least in your specific use case. For what it is worth, the error you were getting was because you need to pass a list to the compositor object...but the rest of your code wouldn't have work anyway. So the easier solution should be...

from satpy.writers import to_image
from satpy.composites import GenericCompositor

r_band = to_image(scn["1"])
r_band.stretch_hist_equalize(True)

g_band = to_image(scn["2"])
g_band.stretch_hist_equalize(True)

b_band = to_image(scn["3"])
b_band.stretch_hist_equalize(True)

comp = GenericCompositor("test1")
test1 = comp((r_band.data, g_band.data，b_band.data))
img = to_image(test1)
img.save("test1.tif")

Note that to_image is returning an XRImage (from the trollimage library). The GenericCompositor needs to be passed xarray DataArray objects so we use the .data propery on each of those XRImage objects.

There are slightly different ways to do this involving YAML but this should work...I think.

[haiyangdaozhang]

Thank you! It works!

But when I want to use it to generate true_color without histogram enhancement, the image is completely white.
I checked r_band.data and found that for the synthesized channel, the maximum value of the data is 100. Should the normal value be 0-1?
So I use get_enhanced_data like this, the image is correct. May I use it correctly?

r_band = get_enhanced_image(scn["1"])
r_band.stretch_hist_equalize(True)

g_band =get_enhanced_image(scn["2"])
g_band.stretch_hist_equalize(True)

b_band = get_enhanced_image(scn["3"])
b_band.stretch_hist_equalize(True)

comp = GenericCompositor("test1")
test1 = comp((r_band.data, g_band.data,b_band.data))
img = to_image(test1)
img.save("test1.tif")

And I tried to use modifiers:

scn1.load(["1", "2"], modifiers=('sunz_corrected',))
scn1.keys()

It shows:

[DataID(name='1', wavelength=WavelengthRange(min=0.62, central=0.645, max=0.67, unit='µm'), resolution=250, calibration=<calibration.reflectance>, modifiers=()),
 DataID(name='2', wavelength=WavelengthRange(min=0.841, central=0.8585, max=0.876, unit='µm'), resolution=250, calibration=<calibration.reflectance>, modifiers=())]

modifiers=() It seems that modifiers don't work.

[djhoese]

One step at a time:

1. Don't use get_enhanced_image. You are then using the configured enhancement for those channels which in this case is not a linear enhancement (it is gamma of 1.5). Instead go back to using to_image, but before calling the hist method, call r_band.crude_stretch(0, 100) or whatever reflectance limit you want to use. The reflectances are typically in the range of ~0% to ~120%. This crude stretch method will linearly normalize the data between 0 and 1 which is what is expected in the img.save method.
2. You likely need to resample your Scene before using the modified data. So new_scn = scn.resample(resampler="native") and then use new_scn when you create the individual RGB channel images.

[haiyangdaozhang]

Thank you for your reminder.

When I use modified data, if I load modified data and normal data separately, two channels with the same name will appear. How can I choose the modified data to create RGB images?

such as :

dq1 = DataQuery(name="1", modifiers=("sunz_corrected", "rayleigh_corrected_crefl"))
dq4 = DataQuery(name="4", modifiers=("sunz_corrected", "rayleigh_corrected_crefl"))
dq3 = DataQuery(name="3", modifiers=("sunz_corrected", "rayleigh_corrected_crefl"))
scn.load([dq1,dq3,dq4,"1"])
scn1 = scn.resample(resampler = "native")
scn1.keys()

It shows:

[DataID(name='1', wavelength=WavelengthRange(min=0.62, central=0.645, max=0.67, unit='µm'), resolution=250, calibration=<calibration.reflectance>, modifiers=()),
 DataID(name='1', wavelength=WavelengthRange(min=0.62, central=0.645, max=0.67, unit='µm'), resolution=250, calibration=<calibration.reflectance>, modifiers=('sunz_corrected', 'rayleigh_corrected_crefl')),
 DataID(name='3', wavelength=WavelengthRange(min=0.459, central=0.469, max=0.479, unit='µm'), resolution=500, calibration=<calibration.reflectance>, modifiers=('sunz_corrected', 'rayleigh_corrected_crefl')),
 DataID(name='4', wavelength=WavelengthRange(min=0.545, central=0.555, max=0.565, unit='µm'), resolution=500, calibration=<calibration.reflectance>, modifiers=('sunz_corrected', 'rayleigh_corrected_crefl'))]

How to select the corrected channel 1 or the default channel 1？

Thank you again.

[djhoese]

You can use that same DataQuery object to do scn1[dq1]. If you want the unmodified then you can just use scn1["1"] which defaults to choosing the least modified version of a product. Or you could make another DataQuery object with dq1_unmod = DataQuery(name="1", modifiers=()) and then scn1[dq1_unmod].

[haiyangdaozhang]

Thank you for your help these days! ^_^

[djhoese]

No problem. I think your original issue has been solved, but if you have further questions feel free to comment here or make a new issue.