pytroll / satpy

Python package for earth-observing satellite data processing
http://satpy.readthedocs.org/en/latest/
GNU General Public License v3.0
1.07k stars 295 forks source link

How to colorize images #459

Closed hany closed 5 years ago

hany commented 6 years ago

Hi there!

Just stumbled on this project (and its dependencies) and getting to know its capabilities. I'm looking for a way to produce images from the GOES-16 ABI that are colorized similar to the outputs here:

https://cdn.star.nesdis.noaa.gov/GOES16/ABI/FD/10/20182861845_GOES16-ABI-FD-10-1808x1808.jpg

Looking at the documentation, it appears adding colormaps are possible, however I'm not sure how to use it when saving datasets with scn.save_datasets(). Any help would be appreciated!

djhoese commented 6 years ago

Hi @hany thanks for the issue. It should be possible, but may not be the easiest from python right now. We are also lacking a lot of documentation for enhancements right now as things are still changing relatively quickly.

You can make your own colormap that matches the trollimage Colormap interface (value -> (R, G, B)): https://trollimage.readthedocs.io/en/latest/colormap.html#trollimage.colormap.Colormap

You could then make your own custom enhancement configuration (called enhancements/generic.yaml or enhancements/abi.yaml) and use the existing enhancements as examples: https://github.com/pytroll/satpy/blob/master/satpy/etc/enhancements/generic.yaml

You would make an entry for the specific band(s) you want to add the colormap for specifying name: C01 or whatever band you were processing. Similar to the true_color_crefl enhancement here: https://github.com/pytroll/satpy/blob/master/satpy/etc/enhancements/generic.yaml#L14 Instead of pointing to the crefl_scaling method you would point to colorize and you can create a literal colormap by setting a palettes argument (kind of like but not exactly like https://github.com/pytroll/satpy/blob/master/satpy/etc/enhancements/generic.yaml#L526). One of my student employees is doing something very similar but it looks like she hasn't pushed her changed to her fork yet. @katherinekolman maybe when you push your branch you could point to your enhancement.yaml file as an example for @hany.

The save_datasets method will automatically pick up on your custom config if you have it in your local directory or point to the root directory by setting a PPP_CONFIG_DIR environment variable. We hope to make this simpler in the future. Hope this helps as a starting point at least.

hany commented 6 years ago

@djhoese thanks very much for your help and explanation here. I created a test colormap and an associated enhancements/abi.yaml file with the following:

enhancements:
  color:
    name: C13
    standard_name: C13
    operations:
    - name: colorize
      method: !!python/name:satpy.enhancements.colorize
      kwargs:
        palettes:
          - {filename: /tmp/binary_colormap.npy, min_value: 223.15, max_value: 303.15}

Taken from the example provided here: https://satpy.readthedocs.io/en/latest/writers.html#colorizing-and-palettizing-using-user-supplied-colormaps

However, when I attempt to save_datasets, it doesn't appear to load my enhancement. Really simply, I'm doing:

BASE_DIR = '.'
all_filenames = [glob(fn.replace('C01', 'C0[123]*')[:len(BASE_DIR) + 50] + '*.nc') for fn in sorted(glob(os.path.join(BASE_DIR, 'OR*-RadF-*C*.nc')))]

scenes = [Scene(reader='abi_l1b', filenames=filenames) for filenames in all_filenames]
print("Number of Scenes: ", len(scenes))

scene = scenes[0]
scene.load(['C13'])
new_scn = scene.resample(scene.min_area(), resampler='native')

with ProgressBar():
    new_scn.save_datasets()

When running with debug, I can see the default enhancement is being picked up:

DEBUG:satpy.writers:Enhancement configuration options: [{'name': 'stretch', 'method': <function stretch at 0x118dec598>, 'kwargs': {'stretch': 'linear'}}]

I went as far as editing the main https://github.com/pytroll/satpy/blob/master/satpy/etc/enhancements/generic.yaml file on disk to add this enhancement, but it's not being picked up. I think I'm missing a step somewhere. Any clues?

djhoese commented 6 years ago

Remove the standard_name part of your config or change it to toa_brightness_temperature. You have it set to C13 which is not what the standard_name for that channel is.

Edit: Note that if you removed the first name: part then the enhancement would be used for all datasets matching standard_name. If you remove standard_name it would match for only datasets matching the exact name. If you have both it must match both of those items. And since you are putting this in the abi.yaml config it will only be used when the sensor of the dataset being enhanced is abi. You could accomplish the same thing in generic.yaml by specifying sensor: abi in addition to name and/or standard_name.

hany commented 6 years ago

Thanks so much @djhoese! That was indeed it. I'm going to continue with this to see if I can get any colored images and report back for any future folk going through the same thing.

hany commented 6 years ago

OK! Finally figured this out and got it working. For anyone else stumbling across this issue:

Thanks for your help!

djhoese commented 6 years ago

@pnuu is this true? The builtin trollimage colormaps don't work?

pnuu commented 6 years ago

Works for me. As an example:


&colorizefun !!python/name:satpy.enhancements.colorize ''
enhancements:
  colorized_ir_clouds:
    standard_name: colorized_ir_clouds
    operations:
      - name: colorize
        method: *colorizefun
        kwargs:
          palettes:
            - {colors: spectral, min_value: 203.15, max_value: 243.149999}
            - {colors: greys, min_value: 243.15, max_value: 303.15}
zxdawn commented 4 years ago

@djhoese @pnuu I met the same problem when colorizing image using enhancement. Here's the simple example of plotting AGRI C12:

import os, glob
from satpy.scene import Scene
from satpy.utils import debug_on

debug_on()
os.environ['PPP_CONFIG_DIR'] = '/yin_raid/xin/satpy_config/'

date = '20190725'
save_dir = './figures/satpy/'
save_path = save_dir + date + '/'

os.makedirs(save_path, exist_ok = True)

filenames = glob.glob('/xin/data/FY4A/20190807/FY4A-_AGRI*4000M_V0001.HDF')

scn = Scene(filenames, reader='agri_l1')

channel = 'C12'
scn.load([channel], generate=False, calibration='brightness_temperature')
img = scn.resample('LEKIMA_1km', cache_dir='./')

img.save_datasets(base_dir = save_path, filename='{sensor}_{name}_{start_time:%Y-%m-%d_%H-%M}.png',
                     compute = True,
                     datasets = [channel],
                     writer = 'simple_image'
                 )

I created the agri.yaml under /yin_raid/xin/satpy_config/enhancements/:

&colorizefun !!python/name:satpy.enhancements.colorize ''

enhancements:
  AGRI_C12:
    sensor: AGRI
    name: C12
    standard_name: toa_brightness_temperature
    operations:
      - name: colorize
        method: *colorizefun
        kwargs:
          palettes:
            - {colors: spectral, min_value: 193.15, max_value: 253.149999}
            - {colors: greys, min_value: 253.15, max_value: 303.15}

Even sensor and standard_name are deleted from the enhancement, save_datasets still uses stretch as the enhancement.

However, if I copy that enhancement to generic.yaml, it works well.

djhoese commented 4 years ago

Can you print out the scn['C12'] and paste the output here? I'm most concerned with the Attributes. Is the sensor exactly agri? I think the sensor in the enhancement config should say agri lowercase.

zxdawn commented 4 years ago

Ah, it's uppercase:

Attributes:
    orbital_parameters:   {'satellite_nominal_latitude': array(0., dtype=floa...
    center_wavelength:    10.8um
    lut_key:              CALChannel12
    resolution:           4000
    long_name:            10.8um channel 4km image data layer
    standard_name:        toa_brightness_temperature
    file_key:             NOMChannel12
    file_type:            agri_l1_4000m
    level:                None
    calibration:          brightness_temperature
    platform_name:        FY-4A
    wavelength:           (10.3, 10.8, 11.1)
    name:                 C12
    sensor:               AGRI
    fill_value:           65535
    units:                K
    polarization:         None
    band_names:           band12(band number is range from 1 to 14)
    valid_range:          [100. 500.]
    modifiers:            ()
    start_time:           2019-08-07 06:00:04.568000
    end_time:             2019-08-07 06:12:43.188000
    area:                 Area ID: DISK\nDescription: AGRI DISK area\nProject...
    ancillary_variables:  []

Should I fix it in the reader?

djhoese commented 4 years ago

I'm not sure if we've standardized it necessarily, but we probably lean towards lowercase for sensors. If you change it in the reader then it should work.

zxdawn commented 4 years ago

Thanks. It works after .lower() is applied to the sensor name. I will push it to that reader: https://github.com/pytroll/satpy/blob/24d5bdbd94d3665bab2ccac8c008efef5ff06acd/satpy/readers/agri_l1.py#L106