From 738348eb2ddbf01b6b90b11fe45932d39256c541 Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Thu, 15 Sep 2016 10:18:46 +0200 Subject: [PATCH] Add imagery capability for OCA cloud parameters. Only the cloud top pressure (ctp) parameters is okay so far. Need to check effective radius and COT --- mpop/imageo/palettes.py | 254 ++++++++++++++++++++++++++++--------- mpop/instruments/seviri.py | 72 ++++++++++- 2 files changed, 266 insertions(+), 60 deletions(-) diff --git a/mpop/imageo/palettes.py b/mpop/imageo/palettes.py index 3205299e..867045ca 100644 --- a/mpop/imageo/palettes.py +++ b/mpop/imageo/palettes.py @@ -1,6 +1,6 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2009, 2013, 2015. +# Copyright (c) 2009, 2013, 2015, 2016. # SMHI, # Folkborgsvägen 1, @@ -30,16 +30,18 @@ """Palette holder module. """ +import numpy as np + def tv_legend(): """Palette for TV. """ legend = [] - legend.append(( 0, 0, 0)) # Unprocessed: Black - legend.append(( 0, 120, 0)) # Land - legend.append(( 0, 0, 215)) # Sea: Blue - legend.append(( 0, 120, 0)) # Land (Snow on land) - legend.append(( 0, 0, 215)) # Sea: Blue (Snow/Ice on sea) + legend.append((0, 0, 0)) # Unprocessed: Black + legend.append((0, 120, 0)) # Land + legend.append((0, 0, 215)) # Sea: Blue + legend.append((0, 120, 0)) # Land (Snow on land) + legend.append((0, 0, 215)) # Sea: Blue (Snow/Ice on sea) for i in range(5, 256): # All other pixel values are grey according to IR temp. @@ -52,9 +54,9 @@ def vv_legend(): """Palette for Swedish road authorities (Vägverket). """ legend = [] - legend.append(( 0, 0, 0)) # Unprocessed: Black - legend.append(( 0, 120, 0)) # Land - legend.append(( 0, 0, 215)) # Sea: Blue + legend.append((0, 0, 0)) # Unprocessed: Black + legend.append((0, 120, 0)) # Land + legend.append((0, 0, 215)) # Sea: Blue # Cloud type values 5 to 8: legend.append((255, 150, 0)) # Very low cumuliform legend.append((255, 100, 0)) # Very low @@ -67,17 +69,19 @@ def vv_legend(): return convert_palette(legend) + def cloud_phase(): """Palette for cloud phase. """ legend = [] - legend.append(( 0, 0, 0)) # Unprocessed: Black - legend.append(( 0, 0, 215)) # Water Clouds: Blue - legend.append(( 240, 240, 240)) # Ice Clouds: Almost White - legend.append(( 120, 120, 0)) # Uncertain Phase: ? - + legend.append((0, 0, 0)) # Unprocessed: Black + legend.append((0, 0, 215)) # Water Clouds: Blue + legend.append((240, 240, 240)) # Ice Clouds: Almost White + legend.append((120, 120, 0)) # Uncertain Phase: ? + return convert_palette(legend) + def cms_modified(): """Palette for regular cloud classification. """ @@ -89,8 +93,8 @@ def nwcsaf_cloudtype(): """ legend = [] legend.append((100, 100, 100)) # Unprocessed: Grey - legend.append(( 0, 120, 0)) - legend.append(( 0, 0, 0)) # Sea: Black + legend.append((0, 120, 0)) + legend.append((0, 0, 0)) # Sea: Black legend.append((250, 190, 250)) # Snow legend.append((220, 160, 220)) # Sea-ice @@ -105,12 +109,12 @@ def nwcsaf_cloudtype(): legend.append((255, 255, 255)) # Very high cumuliform legend.append((230, 230, 230)) # Very high - legend.append(( 0, 80, 215)) # Semi-transparent thin - legend.append(( 0, 180, 230)) # Semi-transparent medium - legend.append(( 0, 240, 240)) # Semi-transparent thick - legend.append(( 90, 200, 160)) # Semi-transparent above + legend.append((0, 80, 215)) # Semi-transparent thin + legend.append((0, 180, 230)) # Semi-transparent medium + legend.append((0, 240, 240)) # Semi-transparent thick + legend.append((90, 200, 160)) # Semi-transparent above legend.append((200, 0, 200)) # Broken - legend.append(( 95, 60, 30)) # Undefined: Brown + legend.append((95, 60, 30)) # Undefined: Brown return convert_palette(legend) @@ -119,7 +123,7 @@ def ctth_height(): """CTTH height palette. """ legend = [] - legend.append(( 0, 0, 0)) + legend.append((0, 0, 0)) legend.append((255, 0, 216)) # 0 meters legend.append((126, 0, 43)) legend.append((153, 20, 47)) @@ -131,11 +135,11 @@ def ctth_height(): legend.append((216, 255, 0)) legend.append((178, 255, 0)) legend.append((153, 255, 0)) - legend.append(( 0, 255, 0)) - legend.append(( 0, 140, 48)) - legend.append(( 0, 178, 255)) - legend.append(( 0, 216, 255)) - legend.append(( 0, 255, 255)) + legend.append((0, 255, 0)) + legend.append((0, 140, 48)) + legend.append((0, 178, 255)) + legend.append((0, 216, 255)) + legend.append((0, 255, 255)) legend.append((238, 214, 210)) legend.append((239, 239, 223)) legend.append((255, 255, 255)) # 10,000 meters @@ -178,19 +182,19 @@ def ctth_height_pps(): legend.append((178, 255, 0)) legend.append((153, 255, 0)) legend.append((153, 255, 0)) - legend.append(( 0, 255, 0)) - legend.append(( 0, 255, 0)) - legend.append(( 0, 255, 0)) - legend.append(( 0, 140, 48)) - legend.append(( 0, 140, 48)) - legend.append(( 0, 178, 255)) - legend.append(( 0, 178, 255)) - legend.append(( 0, 178, 255)) - legend.append(( 0, 216, 255)) - legend.append(( 0, 216, 255)) - legend.append(( 0, 255, 255)) - legend.append(( 0, 255, 255)) - legend.append(( 0, 255, 255)) + legend.append((0, 255, 0)) + legend.append((0, 255, 0)) + legend.append((0, 255, 0)) + legend.append((0, 140, 48)) + legend.append((0, 140, 48)) + legend.append((0, 178, 255)) + legend.append((0, 178, 255)) + legend.append((0, 178, 255)) + legend.append((0, 216, 255)) + legend.append((0, 216, 255)) + legend.append((0, 255, 255)) + legend.append((0, 255, 255)) + legend.append((0, 255, 255)) legend.append((238, 214, 210)) legend.append((238, 214, 210)) legend.append((239, 239, 223)) @@ -290,44 +294,180 @@ def convert_palette(palette): i[1] / 255.0, i[2] / 255.0)) return new_palette - + + def convert_palette2colormap(palette): """Convert palette from [0,255] range to [0,1]. - """ + """ from trollimage.colormap import Colormap - j=0 - n_pal = len(palette)-1 - values=[] - colors=[] + j = 0 + n_pal = len(palette) - 1 + values = [] + colors = [] - red = [r for (r, g, b) in palette] + red = [r for (r, g, b) in palette] green = [g for (r, g, b) in palette] - blue = [b for (r, g, b) in palette] + blue = [b for (r, g, b) in palette] - max_pal = max( max(red), max(blue), max(green) ) + max_pal = max(max(red), max(blue), max(green)) if max_pal <= 1.0: # print "palette already normalized" - denom = 1.0 + denom = 1.0 else: # print "palette normalized to 255" denom = 255.0 for i in palette: - values.append( (n_pal-j) / float(n_pal) ) + values.append((n_pal - j) / float(n_pal)) colors.append((i[0] / denom, i[1] / denom, i[2] / denom)) - j=j+1 + j = j + 1 # reverse order to the entries - values=values[::-1] - colors=colors[::-1] + values = values[::-1] + colors = colors[::-1] - #for i in palette: + # for i in palette: # values.append( j / float(n_pal)) # colors.append((i[0] / 255.0, i[1] / 255.0, i[2] / 255.0)) # j=j+1 # attention: # Colormap(values, colors) uses the second input option of Colormap - # values has to be a list (not a tuple) and - # colors has to be the corresponding list of color tuples + # values has to be a list (not a tuple) and + # colors has to be the corresponding list of color tuples return Colormap(values, colors) + + +class LogColors(object): + + """ + Defines colors to use with `logdata2image` + + """ + + def __init__(self, nodata, zeros, over, breaks): + self.nodata = nodata + self.zeros = zeros + self.over = over + self.breaks = breaks + + def palette(self, N=256): + """ + Build a palette for logarithmic data images. + + """ + + max_value = self.breaks[-1][0] + + pal = np.zeros((N, 3), dtype=np.uint8) + + b_last, rgb_last = self.breaks[0] + for b, rgb in self.breaks[1:]: + # Get a slice of the palette array for the current interval + p = pal[ + np.log(b_last + 1) * N / np.log(max_value):np.log(b + 1) * N / np.log(max_value)] + for i in range(3): # red, green, blue + p[:, i] = np.linspace(rgb_last[i], rgb[i], p.shape[0]) + b_last = b + rgb_last = rgb + + pal[0] = self.nodata + pal[1] = self.zeros + pal[-1] = self.over + + return pal + + +class TriColors(LogColors): + + """ + Use three color tones in the intervals between the elements of *breaks*. + + """ + color_tones = [((0, 0, 200), (150, 150, 255)), # dark to light blue + ((150, 150, 0), (255, 255, 8)), # greyish to bright yellow + ((230, 150, 100), (230, 0, 0))] # green to red + + nodata = (0, 0, 0) # black + # zeros = (20, 0, 20) # dark purple + # black #There is no need to mark zeros with another col + zeros = (0, 0, 0) + over = (255, 0, 0) # bright red + + def __init__(self, breaks): + breaks = [(breaks[0], TriColors.color_tones[0][0]), + (breaks[1], TriColors.color_tones[0][1]), + + (breaks[1], TriColors.color_tones[1][0]), + (breaks[2], TriColors.color_tones[1][1]), + + (breaks[2], TriColors.color_tones[2][0]), + (breaks[3], TriColors.color_tones[2][1])] + + LogColors.__init__(self, TriColors.nodata, TriColors.zeros, + TriColors.over, breaks) + +CPP_COLORS = {'cpp_cot': TriColors([0, 3.6, 23, 700]), # ISCCP intervals + 'cpp_reff': TriColors([0, 10, 20, 1000])} + +CPP_COLORS['cot'] = CPP_COLORS['cpp_cot'] +CPP_COLORS['reff'] = CPP_COLORS['cpp_reff'] + + +def get_ctp_legend(): + """ + Get the Cloud Top Pressure color palette + """ + + legend = [] + legend.append((0, 0, 0)) # No data + legend.append((255, 0, 216)) # 0: 1000-1050 hPa (=100000-105000 Pa) + legend.append((126, 0, 43)) # 1: 950-1000 hPa + legend.append((153, 20, 47)) # 2: 900-950 hPa + legend.append((178, 51, 0)) # 3: 850-900 hPa + legend.append((255, 76, 0)) # 4: 800-850 hPa + legend.append((255, 102, 0)) # 5: 750-800 hPa + legend.append((255, 164, 0)) # 6: 700-750 hPa + legend.append((255, 216, 0)) # 7: 650-700 hPa + legend.append((216, 255, 0)) # 8: 600-650 hPa + legend.append((178, 255, 0)) # 9: 550-600 hPa + legend.append((153, 255, 0)) # 10: 500-550 hPa + legend.append((0, 255, 0)) # 11: 450-500 hPa + legend.append((0, 140, 48)) # 12: 400-450 hPa + legend.append((0, 178, 255)) # 13: 350-400 hPa + legend.append((0, 216, 255)) # 14: 300-350 hPa + legend.append((0, 255, 255)) # 15: 250-300 hPa + legend.append((238, 214, 210)) # 16: 200-250 hPa + legend.append((239, 239, 223)) # 17: 150-200 hPa + legend.append((255, 255, 255)) # 18: 100-150 hPa + legend.append((255, 255, 255)) # 19: 50-100 hPa + legend.append((255, 255, 255)) # 20: 0-50 hPa (=0-5000 Pa) + + palette = convert_palette(legend) + return palette + + +def get_reff_legend(): + return get_log_legend('reff') + + +def get_cot_legend(): + return get_log_legend('cot') + + +def get_log_legend(product_name): + # This is the same data as is used in logdata2image (when indata as for + # the calls from cppimage) + return CPP_COLORS[product_name].palette() + + +def oca_get_scenetype_legend(): + + # Colorize using PPS/CPP palette + legend = np.array([[170, 130, 255], # purple/blue for liquid (cph == 1) + [220, 200, 255], # almost white for ice (cph == 2) + [255, 200, 200] # Redish for multi layer clouds + ]) + legend = np.vstack([np.zeros((111, 3)), legend]) + palette = convert_palette(legend) + return palette diff --git a/mpop/instruments/seviri.py b/mpop/instruments/seviri.py index 2f04e53f..c09dcbe4 100644 --- a/mpop/instruments/seviri.py +++ b/mpop/instruments/seviri.py @@ -31,14 +31,49 @@ LOG = logging.getLogger(__name__) -import os.path - try: from pyorbital.astronomy import sun_zenith_angle as sza except ImportError: sza = None +from mpop.imageo import palettes + +oca_palette_func = {'ll_ctp': palettes.get_ctp_legend, + 'ul_ctp': palettes.get_ctp_legend, + 'ul_cot': palettes.get_cot_legend, + 'll_cot': palettes.get_cot_legend, + 'reff': palettes.get_reff_legend, + 'scenetype': palettes.oca_get_scenetype_legend} + + +def _arrange_log_data(arr, max_value, no_data_value): + """ + Prepare logarithmic data for creating an image. + + """ + MAX_IM_VAL = 2**8 - 1 + + # Logarithmic data should never be negative + assert ((arr >= 0) + (arr == no_data_value)).all(), \ + "Negative values encountered in cloud optical thickness" + + # Confine image data values to the range [2, MAX_IM_VAL] + arr_log = np.log(arr + 1.) # arr == 0 -> arr_log = 0 + cot_im_data = arr_log * (MAX_IM_VAL - 3) / np.log(max_value + 1.) + 2. + cot_im_data[cot_im_data > MAX_IM_VAL] = MAX_IM_VAL + + # Now that the data is adjusted, cast it to uint8 ([0, 2**8 - 1]) + cot_im_data = cot_im_data.astype(np.uint8) + + # Give no-data values a special color + cot_im_data[arr == no_data_value] = 0 + # Give arr == 0 a special color + cot_im_data[arr == 0] = 1 + + return cot_im_data + + class SeviriCompositer(VisirCompositer): """This class sets up the Seviri instrument channel list. @@ -337,7 +372,7 @@ def day_microphysics(self, wintertime=False, fill_value=None): self.area, self.time_slot, crange=crange, - fill_value=fill_value, mode="RGB") + fill_value=fill_value, mode="RGB") if wintertime: img.gamma((1.0, 1.5, 1.0)) else: @@ -347,3 +382,34 @@ def day_microphysics(self, wintertime=False, fill_value=None): day_microphysics.prerequisites = refl39_chan.prerequisites | set( [0.8, 10.8]) + + def oca(self, fieldname): + """Make an OCA cloud parameter image""" + + palette = oca_palette_func[fieldname]() + data = getattr(getattr(self['OCA'], fieldname), 'data') + if fieldname in ['ul_ctp', 'll_ctp']: + data = (22. - data / 5000.).astype('Int16') + + elif fieldname in ['reff']: + data = (data * 1000000. + 0.5).astype('uint8') + + elif fieldname in ['ul_cot', 'll_cot']: + import pdb + max_value = palettes.CPP_COLORS['cot'].breaks[-1][0] + no_data = -1000 # FIXME! + data = _arrange_log_data(data, max_value, no_data) + pdb.set_trace() + + else: + raise NotImplementedError( + "No imagery for parameter %s implemented yet...", fieldname) + + img = geo_image.GeoImage(data, self.area, + self.time_slot, + fill_value=(0), + mode="P", + palette=palette) + return img + + oca.prerequisites = set(['OCA'])