python-mpop 1.0.0-1 / usr / lib / python2.7 / dist-packages / mpop / scene.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2010, 2011, 2012, 2013.

# Author(s):
 
#   Martin Raspaud <martin.raspaud@smhi.se>
#   Esben S. Nielsen <esn@dmi.dk>

# This file is part of mpop.

# mpop is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.

# mpop is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.

# You should have received a copy of the GNU General Public License along with
# mpop.  If not, see <http://www.gnu.org/licenses/>.

"""The :mod:`mpop.scene` module defines satellite scenes. They are defined as
generic classes, to be inherited when needed.

A scene is a set of :mod:`mpop.channel` objects for a given time, and sometimes
also for a given area.
"""
import ConfigParser
import copy
import datetime
import os.path
import types
import weakref
import sys

import numpy as np

from mpop import CONFIG_PATH
from mpop.channel import Channel, NotLoadedError
from mpop.logger import LOG
from mpop.utils import OrderedConfigParser


try:
    # Work around for on demand import of pyresample. pyresample depends 
    # on scipy.spatial which memory leaks on multiple imports
    is_pyresample_loaded = False
    from pyresample.geometry import AreaDefinition, SwathDefinition
    import mpop.projector
    is_pyresample_loaded = True
except ImportError:
    LOG.warning("pyresample missing. Can only work in satellite projection")
    

class Satellite(object):
    """This is the satellite class. It contains information on the satellite.
    """

    def __init__(self, (satname, number, variant)=(None, None, None)):
        try:
            self.satname = satname or "" or self.satname
        except AttributeError:
            self.satname = satname or ""
        try:
            self.number = number or "" or self.number
        except AttributeError:
            self.number = number or ""
        try:
            self.variant = variant or "" or self.variant
        except AttributeError:
            self.variant = variant or ""

    @property
    def fullname(self):
        """Full name of the satellite, that is platform name and number
        (eg "metop02").
        """
        return self.variant + self.satname + self.number

    @classmethod
    def remove_attribute(cls, name):
        """Remove an attribute from the class.
        """
        return delattr(cls, name)

    @classmethod
    def add_method(cls, func):
        """Add a method to the class.
        """
        return setattr(cls, func.__name__, func)

    def add_method_to_instance(self, func):
        """Add a method to the instance.
        """
        return setattr(self, func.__name__,
                       types.MethodType(func, self.__class__))


class SatelliteScene(Satellite):
    """This is the satellite scene class. It is a capture of the satellite
    (channels) data at given *time_slot* and *area_id*/*area*.
    """

    def __init__(self, time_slot=None, area_id=None, area=None,
                 orbit=None, satellite=(None, None, None)):

        Satellite.__init__(self, satellite)
        
        if(time_slot is not None and
           not isinstance(time_slot, datetime.datetime)):
            raise TypeError("Time_slot must be a datetime.datetime instance.")
        
        self.time_slot = time_slot


        self.area_id = None
        self.area_def = None
        
        if(area_id is not None):
            from warnings import warn
            warn("The *area_id* attribute is deprecated."
                 "Please use *area* instead.",
                 DeprecationWarning)
            if(not isinstance(area_id, str)):
                raise TypeError("Area must be a string.")

        self.area = area_id

        if area is not None:
            self.area = area
        
        if(orbit is not None and
           not isinstance(orbit, str)):
            raise TypeError("Orbit must be a string.")
        
        self.orbit = orbit


        self.info = {}
        
        self.lat = None
        self.lon = None


    def get_area(self):
        """Getter for area.
        """
        return self.area_def or self.area_id

    def set_area(self, area):
        """Setter for area.
        """
        if (area is None):
            self.area_def = None
            self.area_id = None
        elif(isinstance(area, str)):
            self.area_id = area
            self.area_def = None
        else:
            try:
                dummy = area.area_extent
                dummy = area.x_size
                dummy = area.y_size
                dummy = area.proj_id
                dummy = area.proj_dict
                self.area_def = area
                self.area_id = None
            except AttributeError:
                try:
                    dummy = area.lons
                    dummy = area.lats
                    self.area_def = area
                    self.area_id = None
                except AttributeError:
                    raise TypeError(("Malformed area argument. "
                                    "Should be a string or an area object. "
                                    "Not %s") % type(area))

    area = property(get_area, set_area)

class SatelliteInstrumentScene(SatelliteScene):
    """This is the satellite instrument class. It is an abstract channel
    container, from which all concrete satellite scenes should be derived.

    The constructor accepts as optional arguments the *time_slot* of the scene,
    the *area* on which the scene is defined (this can be use for slicing of
    big datasets, or can be set automatically when loading), and *orbit* which
    is a string giving the orbit number.
    """
    channel_list = []

    def __init__(self, time_slot=None, area_id=None, area=None,
                 orbit=None, satellite=(None, None, None), instrument=None):

        SatelliteScene.__init__(self, time_slot, area_id, area,
                                orbit, satellite)
        
        try:
            self.instrument_name = instrument or self.instrument_name
        except AttributeError:
            self.instrument_name = None
            
        self.channels = []

        try:
            conf = OrderedConfigParser()
            conf.read(os.path.join(CONFIG_PATH, self.fullname+".cfg"))

            for section in conf.sections():
                if(not section[:-1].endswith("level") and
                   not section.endswith("granules") and
                   section.startswith(self.instrument_name)):
                    name = eval(conf.get(section, "name"))
                    try:
                        w_range = eval(conf.get(section, "frequency"))
                    except ConfigParser.NoOptionError:
                        w_range = (-np.inf, -np.inf, -np.inf)
                    try:
                        resolution = eval(conf.get(section, "resolution"))
                    except ConfigParser.NoOptionError:
                        resolution = 0
                    self.channels.append(Channel(name=name,
                                                 wavelength_range=w_range,
                                                 resolution=resolution))

        except (ConfigParser.NoSectionError, ConfigParser.NoOptionError):
            for name, w_range, resolution in self.channel_list:
                self.channels.append(Channel(name=name,
                                             wavelength_range=w_range,
                                             resolution=resolution))

        self.channels_to_load = set([])

    def __getitem__(self, key, aslist = False):
        if(isinstance(key, float)):
            channels = [chn for chn in self.channels
                        if(hasattr(chn, "wavelength_range") and
                           chn.wavelength_range[0] <= key and
                           chn.wavelength_range[2] >= key)]
            channels = sorted(channels,
                              lambda ch1,ch2:
                                  ch1.__cmp__(ch2, key))
            
        elif(isinstance(key, str)):
            channels = [chn for chn in self.channels
                        if chn.name == key]
            channels = sorted(channels)

        elif(isinstance(key, int)):
            channels = [chn for chn in self.channels
                        if int(np.round(chn.resolution)) == key]
            channels = sorted(channels)

        elif(isinstance(key, (tuple, list))):
            if len(key) == 0:
                raise KeyError("Key list must contain at least one element.")
            channels = self.__getitem__(key[0], aslist = True)
            if(len(key) > 1 and len(channels) > 0):
                dummy_instance = SatelliteInstrumentScene()
                dummy_instance.channels = channels
                channels = dummy_instance.__getitem__(key[1:], aslist = True)
        else:
            raise TypeError("Malformed key: " + str(key))

        if len(channels) == 0:
            raise KeyError("No channel corresponding to "+str(key)+".")
        elif aslist:
            return channels
        else:
            return channels[0]

    def __setitem__(self, key, data):
        # Add a channel if it is not already in the scene. Works only if key is
        # a string.
        try:
            if key not in self:
                # if it's a blob with name and data, add it as is.
                if hasattr(data, "name") and hasattr(data, "data"):
                    self.channels.append(data)
                else:
                    kwargs = {"name": key}
                    for attr in ["wavelength_range", "resolution"]:
                        try:
                            kwargs[attr] = getattr(data, attr)
                        except (AttributeError, NameError):
                            pass
                    self.channels.append(Channel(**kwargs))
        except AttributeError:
            pass

        # Add the data.
        if isinstance(data, np.ma.core.MaskedArray):
            self[key].data = data
        else:
            try:
                self[key].data = data.data
            except AttributeError:
                self[key].data = data



        # if isinstance(data, Channel):
        #     self.channels.append(Channel(name=key,
        #                              wavelength_range=data.wavelength_range,
        #                              resolution=data.resolution))
        #     self[key].data = data.data
        # else:
        #     try:
        #         self[key].data = data
        #     except KeyError:
        #         self.channels.append(Channel(name=key))
        #         self[key].data = data                             
                


    def __str__(self):
        return "\n".join([str(chn) for chn in self.channels])

    def __iter__(self):
        return self.channels.__iter__()


    def _set_reader(self, pformat):
        """Gets the reader for *pformat* format, and puts it in the `reader`
        attribute.
        """
       
        elements = pformat.split(".")
        if len(elements) == 1:
            reader_module = pformat
            reader = "mpop.satin."+reader_module

            # Loading old style plugins
            reader_module = pformat
            LOG.info("old style plugin: " + pformat)
            try:
                # Look for builtin reader
                loader = __import__(reader, globals(),
                                    locals(), ['load'])
            except ImportError:
                # Look for custom reader
                loader = __import__(reader_module, globals(),
                                    locals(), ['load'])

            # Build a custom Reader plugin on the fly...
            from mpop.plugin_base import Reader
            reader_class = type(elements[-1].capitalize() + "Reader",
                                (Reader,),
                                {"pformat": elements[-1]})

            reader_instance = reader_class(self)

            # ... and set its "load" attribute with the "load" function of the
            # loader module
            loader = getattr(loader, "load")
            setattr(reader_instance, "load", loader)

            setattr(self, elements[-1] + "_reader", reader_instance)


        else:
            reader_module = ".".join(elements[:-1])
            reader_class = elements[-1]
        
            reader = "mpop.satin."+reader_module
            try:
                # Look for builtin reader
                loader = __import__(reader, globals(),
                                    locals(), [reader_class])
            except ImportError:
                # Look for custom reader
                loader = __import__(reader_module, globals(),
                                    locals(), [reader_class])
            loader = getattr(loader, reader_class)
            reader_instance = loader(self)
            setattr(self, loader.pformat + "_reader", reader_instance)

        return reader_instance





    def load(self, channels=None, load_again=False, area_extent=None, **kwargs):
        """Load instrument data into the *channels*. *Channels* is a list or a
        tuple containing channels we will load data into, designated by there
        center wavelength (float), resolution (integer) or name (string). If
        None, all channels are loaded.

        The *load_again* boolean flag allows to reload the channels even they
        have already been loaded, to mirror changes on disk for example. This
        is false by default.

        The *area_extent* keyword lets you specify which part of the data to
        load. Given as a 4-element sequence, it defines the area extent to load
        in satellite projection.

        The other keyword arguments are passed as is to the reader
        plugin. Check the corresponding documentation for more details.
        """

        # Set up the list of channels to load.
        if channels is None:
            for chn in self.channels:
                self.channels_to_load |= set([chn.name])

        elif(isinstance(channels, (list, tuple, set))):
            self.channels_to_load = set()
            for chn in channels:
                try:
                    self.channels_to_load |= set([self[chn].name])
                except KeyError:
                    self.channels_to_load |= set([chn])

        else:
            raise TypeError("Channels must be a list/"
                            "tuple/set of channel keys!")

        loaded_channels = [chn.name for chn in self.loaded_channels()]
        if load_again:
            for chn in self.channels_to_load:
                if chn in loaded_channels:
                    self.unload(chn)
                    loaded_channels = []
        else:
            for chn in loaded_channels:
                self.channels_to_load -= set([chn])

        # find the plugin to use from the config file
        conf = ConfigParser.ConfigParser()
        try:
            conf.read(os.path.join(CONFIG_PATH, self.fullname + ".cfg"))
            if len(conf.sections()) == 0:
                raise ConfigParser.NoSectionError(("Config file did "
                                                    "not make sense"))
            levels = [section for section in conf.sections()
                      if section.startswith(self.instrument_name+"-level")]
        except ConfigParser.NoSectionError:
            LOG.warning("Can't load data, no config file for " + self.fullname)
            self.channels_to_load = set()
            return
        
        levels.sort()

        if levels[0] == self.instrument_name+"-level1":
            levels = levels[1:]

        if len(levels) == 0:
            raise ConfigParser.NoSectionError(
                self.instrument_name+"-levelN (N>1) to tell me how to"+
                " read data... Not reading anything.")

        for level in levels:
            if len(self.channels_to_load) == 0:
                return

            LOG.debug("Looking for sources in section "+level)
            reader_name = conf.get(level, 'format')
            try:
                reader_name = eval(reader_name)
            except NameError:
                reader_name = str(reader_name)
            LOG.debug("Using plugin mpop.satin."+reader_name)

            # read the data
            reader = "mpop.satin."+reader_name
            try:
                reader_instance = self._set_reader(reader_name)
                if area_extent is not None:
                    if(isinstance(area_extent, (tuple, list)) and
                       len(area_extent) == 4):
                        kwargs["area_extent"] = area_extent
                    else:
                        raise ValueError("Area extent must be a sequence of "
                                         "four numbers.")
                    
                reader_instance.load(self, **kwargs)
            except ImportError, e:
                LOG.exception("ImportError while loading "+reader_name+": "
                              + str(e))
                continue
            loaded_channels = set([chn.name for chn in self.loaded_channels()])
            just_loaded = loaded_channels & self.channels_to_load
            if len(just_loaded) == 0:
                LOG.info("No channels loaded with " + reader + ".")
            self.channels_to_load -= loaded_channels
            LOG.debug("Successfully loaded: "+str(just_loaded))
            
        if len(self.channels_to_load) > 0:
            LOG.warning("Unable to import channels "
                        + str(self.channels_to_load))

        self.channels_to_load = set()

    def save(self, filename, to_format="netcdf4", **options):
        """Saves the current scene into a file of format *to_format*. Supported
        formats are:
        
        - *netcdf4*: NetCDF4 with CF conventions.
        """

        writer = "satout." + to_format
        try:
            writer_module = __import__(writer, globals(), locals(), ["save"])
        except ImportError, err:
            raise ImportError("Cannot load "+writer+" writer: "+str(err))

        return writer_module.save(self, filename, **options)

    def unload(self, *channels):
        """Unloads *channels* from
        memory. :meth:`mpop.scene.SatelliteInstrumentScene.load` must be called
        again to reload the data.
        """
        for chn in channels:
            try:
                self[chn].data = None
            except AttributeError:
                LOG.warning("Can't unload channel" + str(chn))
        
        
    def add_to_history(self, message):
        """Adds a message to history info.
        """
        timestr = datetime.datetime.utcnow().isoformat()
        timed_message = str(timestr + " - " + message)
        if not self.info.get("history", ""):
            self.info["history"] = timed_message
        else:
            self.info["history"] += "\n" + timed_message
            

    def check_channels(self, *channels):
        """Check if the *channels* are loaded, raise an error otherwise.
        """
        for chan in channels:
            if not self[chan].is_loaded():
                raise NotLoadedError("Required channel %s not loaded,"
                                     " aborting."%chan)

        return True

    def loaded_channels(self):
        """Return the set of loaded_channels.
        """
        return set([chan for chan in self.channels if chan.is_loaded()])

    def project(self, dest_area, channels=None, precompute=False, mode=None, radius=None):
        """Make a copy of the current snapshot projected onto the
        *dest_area*. Available areas are defined in the region configuration
        file (ACPG). *channels* tells which channels are to be projected, and
        if None, all channels are projected and copied over to the return
        snapshot.

        If *precompute* is set to true, the projecting data is saved on disk
        for reusage. *mode* sets the mode to project in: 'quick' which works
        between cartographic projections, and, as its denomination indicates,
        is quick (but lower quality), and 'nearest' which uses nearest
        neighbour for best projection. A *mode* set to None uses 'quick' when
        possible, 'nearest' otherwise.

        *radius* defines the radius of influence for neighbour search in
         'nearest' mode. Setting it to None, or omitting it will fallback to
         default values (5 times the channel resolution) or 10km if the
         resolution is not available.

        Note: channels have to be loaded to be projected, otherwise an
        exception is raised.
        """
        
        if not is_pyresample_loaded:
            # Not much point in proceeding then 
            return self
        
        _channels = set([])

        if channels is None:
            for chn in self.loaded_channels():
                _channels |= set([chn])

        elif(isinstance(channels, (list, tuple, set))):
            for chn in channels:
                try:
                    _channels |= set([self[chn]])
                except KeyError:
                    LOG.warning("Channel "+str(chn)+" not found,"
                                "thus not projected.")
        else:
            raise TypeError("Channels must be a list/"
                            "tuple/set of channel keys!")


        res = copy.copy(self)

        if isinstance(dest_area, str):
            dest_area = mpop.projector.get_area_def(dest_area)

        
        res.area = dest_area
        res.channels = []

        if not _channels <= self.loaded_channels():
            LOG.warning("Cannot project nonloaded channels: %s."
                        %(_channels - self.loaded_channels()))
            LOG.info("Will project the other channels though.")
            _channels = _channels and self.loaded_channels()
        
        cov = {}

        for chn in _channels:
            if chn.area is None:
                if self.area is None:
                    area_name = ("swath_" + self.fullname + "_" +
                                 str(self.time_slot) + "_"
                                 + str(chn.shape))
                    chn.area = area_name
                else:
                    if is_pyresample_loaded:
                        try:                            
                            chn.area = AreaDefinition(
                                self.area.area_id + str(chn.shape),
                                self.area.name,
                                self.area.proj_id,
                                self.area.proj_dict,
                                chn.shape[1],
                                chn.shape[0],
                                self.area.area_extent,
                                self.area.nprocs)
    
                        except AttributeError:
                            try:
                                dummy = self.area.lons
                                dummy = self.area.lats
                                chn.area = self.area
                                area_name = ("swath_" + self.fullname + "_" +
                                             str(self.time_slot) + "_"
                                             + str(chn.shape))
                                chn.area.area_id = area_name
                            except AttributeError:
                                chn.area = self.area + str(chn.shape)
                    else:
                        chn.area = self.area + str(chn.shape)
            else: #chn.area is not None
                if (is_pyresample_loaded and
                    isinstance(chn.area, SwathDefinition) and
                    (not hasattr(chn.area, "area_id") or
                     not chn.area.area_id)):
                    area_name = ("swath_" + self.fullname + "_" +
                                 str(self.time_slot) + "_"
                                 + str(chn.shape) + "_"
                                 + str(chn.name))
                    chn.area.area_id = area_name

            if chn.area == dest_area:
                res.channels.append(chn)
            else:
                if isinstance(chn.area, str):
                    area_id = chn.area
                else:
                    area_id = chn.area_id or chn.area.area_id
                
                if area_id not in cov:
                    if radius is None:
                        if chn.resolution > 0:
                            radius = 5 * chn.resolution
                        else:
                            radius = 10000
                    cov[area_id] = mpop.projector.Projector(chn.area,
                                                            dest_area,
                                                            mode=mode,
                                                            radius=radius)
                    if precompute:
                        try:
                            cov[area_id].save()
                        except IOError:
                            LOG.exception("Could not save projection.")

                try:
                    res.channels.append(chn.project(cov[area_id]))
                except NotLoadedError:
                    LOG.warning("Channel "+str(chn.name)+" not loaded, "
                                "thus not projected.")
        
        # Compose with image object
        try:
            if res._CompositerClass is not None:
                # Pass weak ref to compositor to allow garbage collection
                res.image = res._CompositerClass(weakref.proxy(res))
        except AttributeError:
            pass
        
        return res

if sys.version_info < (2, 5):
    def any(iterable):
        for element in iterable:
            if element:
                return True
        return False


def assemble_segments(segments):
    """Assemble the scene objects listed in *segment_list* and returns the
    resulting scene object.
    """
    from mpop.satellites import GenericFactory
    
    channels = set([])
    for seg in segments:
        channels |= set([chn.name for chn in seg.loaded_channels()])

    seg = segments[0]
    
    new_scene = GenericFactory.create_scene(seg.satname, seg.number,
                                            seg.instrument_name, seg.time_slot,
                                            seg.orbit, variant=seg.variant)
    
    swath_definitions = {}

    for chn in channels:
        new_scene[chn] = np.ma.concatenate([seg[chn].data
                                            for seg in segments
                                            if seg[chn].is_loaded()])
        try:

            area_names = tuple([seg[chn].area.area_id
                                for seg in segments
                                if seg[chn].is_loaded()])
            if area_names not in swath_definitions:
            
                lons = np.ma.concatenate([seg[chn].area.lons[:]
                                          for seg in segments
                                          if seg[chn].is_loaded()])
                lats = np.ma.concatenate([seg[chn].area.lats[:]
                                          for seg in segments
                                          if seg[chn].is_loaded()])
                new_scene[chn].area = SwathDefinition(lons=lons, lats=lats)
                area_name = "+".join(area_names)
                new_scene[chn].area.area_id = area_name
                new_scene[chn].area_id = area_name
                swath_definitions[area_names] = new_scene[chn].area
            else:
                new_scene[chn].area = swath_definitions[area_names]
                new_scene[chn].area_id = new_scene[chn].area.area_id
        except AttributeError:
            pass

    try:
        lons = np.ma.concatenate([seg.area.lons[:] for seg in segments])
        lats = np.ma.concatenate([seg.area.lats[:] for seg in segments])
        new_scene.area = SwathDefinition(lons=lons, lats=lats)
        for chn in channels:
            if any([seg[chn].area for seg in segments]):
                try:
                    lon_arrays = []
                    lat_arrays = []
                    for seg in segments:
                        if seg[chn].area is not None:
                            lon_arrays.append(seg[chn].area.lons[:])
                            lat_arrays.append(seg[chn].area.lats[:])
                        else:
                            lon_arrays.append(seg.area.lons[:])
                            lat_arrays.append(seg.area.lats[:])
                    lons = np.ma.concatenate(lon_arrays)
                    lats = np.ma.concatenate(lat_arrays)
                    new_scene[chn].area = SwathDefinition(lons=lons, lats=lats)
                except AttributeError:
                    pass
    except AttributeError:
        pass

    return new_scene