Advertisement
Guest User

Untitled

a guest
Aug 12th, 2013
77
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 29.34 KB | None | 0 0
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-
  3. # Copyright (c) 2010, 2011, 2012, 2013.
  4.  
  5. # Author(s):
  6.  
  7. #   Martin Raspaud <martin.raspaud@smhi.se>
  8. #   Esben S. Nielsen <esn@dmi.dk>
  9.  
  10. # This file is part of mpop.
  11.  
  12. # mpop is free software: you can redistribute it and/or modify it under the
  13. # terms of the GNU General Public License as published by the Free Software
  14. # Foundation, either version 3 of the License, or (at your option) any later
  15. # version.
  16.  
  17. # mpop is distributed in the hope that it will be useful, but WITHOUT ANY
  18. # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  19. # A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  20.  
  21. # You should have received a copy of the GNU General Public License along with
  22. # mpop.  If not, see <http://www.gnu.org/licenses/>.
  23.  
  24. """The :mod:`mpop.scene` module defines satellite scenes. They are defined as
  25. generic classes, to be inherited when needed.
  26.  
  27. A scene is a set of :mod:`mpop.channel` objects for a given time, and sometimes
  28. also for a given area.
  29. """
  30. import ConfigParser
  31. import copy
  32. import datetime
  33. import os.path
  34. import types
  35. import weakref
  36. import sys
  37.  
  38. import numpy as np
  39.  
  40. from mpop import CONFIG_PATH
  41. from mpop.channel import Channel, NotLoadedError
  42. from mpop.logger import LOG
  43. from mpop.utils import OrderedConfigParser
  44.  
  45.  
  46. try:
  47.     # Work around for on demand import of pyresample. pyresample depends
  48.     # on scipy.spatial which memory leaks on multiple imports
  49.     is_pyresample_loaded = False
  50.     from pyresample.geometry import AreaDefinition, SwathDefinition
  51.     import mpop.projector
  52.     is_pyresample_loaded = True
  53. except ImportError:
  54.     LOG.warning("pyresample missing. Can only work in satellite projection")
  55.    
  56.  
  57. class Satellite(object):
  58.     """This is the satellite class. It contains information on the satellite.
  59.    """
  60.  
  61.     def __init__(self, (satname, number, variant)=(None, None, None)):
  62.         try:
  63.             self.satname = satname or "" or self.satname
  64.         except AttributeError:
  65.             self.satname = satname or ""
  66.         try:
  67.             self.number = number or "" or self.number
  68.         except AttributeError:
  69.             self.number = number or ""
  70.         try:
  71.             self.variant = variant or "" or self.variant
  72.         except AttributeError:
  73.             self.variant = variant or ""
  74.  
  75.     @property
  76.     def fullname(self):
  77.         """Full name of the satellite, that is platform name and number
  78.        (eg "metop02").
  79.        """
  80.         return self.variant + self.satname + self.number
  81.  
  82.     @classmethod
  83.     def remove_attribute(cls, name):
  84.         """Remove an attribute from the class.
  85.        """
  86.         return delattr(cls, name)
  87.  
  88.     @classmethod
  89.     def add_method(cls, func):
  90.         """Add a method to the class.
  91.        """
  92.         return setattr(cls, func.__name__, func)
  93.  
  94.     def add_method_to_instance(self, func):
  95.         """Add a method to the instance.
  96.        """
  97.         return setattr(self, func.__name__,
  98.                        types.MethodType(func, self.__class__))
  99.  
  100.  
  101. class SatelliteScene(Satellite):
  102.     """This is the satellite scene class. It is a capture of the satellite
  103.    (channels) data at given *time_slot* and *area_id*/*area*.
  104.    """
  105.  
  106.     def __init__(self, time_slot=None, area_id=None, area=None,
  107.                  orbit=None, satellite=(None, None, None)):
  108.  
  109.         Satellite.__init__(self, satellite)
  110.        
  111.         if(time_slot is not None and
  112.            not isinstance(time_slot, datetime.datetime)):
  113.             raise TypeError("Time_slot must be a datetime.datetime instance.")
  114.        
  115.         self.time_slot = time_slot
  116.  
  117.  
  118.         self.area_id = None
  119.         self.area_def = None
  120.        
  121.         if(area_id is not None):
  122.             from warnings import warn
  123.             warn("The *area_id* attribute is deprecated."
  124.                  "Please use *area* instead.",
  125.                  DeprecationWarning)
  126.             if(not isinstance(area_id, str)):
  127.                 raise TypeError("Area must be a string.")
  128.  
  129.         self.area = area_id
  130.  
  131.         if area is not None:
  132.             self.area = area
  133.        
  134.         if(orbit is not None and
  135.            not isinstance(orbit, str)):
  136.             raise TypeError("Orbit must be a string.")
  137.        
  138.         self.orbit = orbit
  139.  
  140.  
  141.         self.info = {}
  142.        
  143.         self.lat = None
  144.         self.lon = None
  145.  
  146.  
  147.     def get_area(self):
  148.         """Getter for area.
  149.        """
  150.         return self.area_def or self.area_id
  151.  
  152.     def set_area(self, area):
  153.         """Setter for area.
  154.        """
  155.         if (area is None):
  156.             self.area_def = None
  157.             self.area_id = None
  158.         elif(isinstance(area, str)):
  159.             self.area_id = area
  160.             self.area_def = None
  161.         else:
  162.             try:
  163.                 dummy = area.area_extent
  164.                 dummy = area.x_size
  165.                 dummy = area.y_size
  166.                 dummy = area.proj_id
  167.                 dummy = area.proj_dict
  168.                 self.area_def = area
  169.                 self.area_id = None
  170.             except AttributeError:
  171.                 try:
  172.                     dummy = area.lons
  173.                     dummy = area.lats
  174.                     self.area_def = area
  175.                     self.area_id = None
  176.                 except AttributeError:
  177.                     raise TypeError(("Malformed area argument. "
  178.                                     "Should be a string or an area object. "
  179.                                     "Not %s") % type(area))
  180.  
  181.     area = property(get_area, set_area)
  182.  
  183. class SatelliteInstrumentScene(SatelliteScene):
  184.     """This is the satellite instrument class. It is an abstract channel
  185.    container, from which all concrete satellite scenes should be derived.
  186.  
  187.    The constructor accepts as optional arguments the *time_slot* of the scene,
  188.    the *area* on which the scene is defined (this can be use for slicing of
  189.    big datasets, or can be set automatically when loading), and *orbit* which
  190.    is a string giving the orbit number.
  191.    """
  192.     channel_list = []
  193.  
  194.     def __init__(self, time_slot=None, area_id=None, area=None,
  195.                  orbit=None, satellite=(None, None, None), instrument=None):
  196.  
  197.         SatelliteScene.__init__(self, time_slot, area_id, area,
  198.                                 orbit, satellite)
  199.        
  200.         try:
  201.             self.instrument_name = instrument or self.instrument_name
  202.         except AttributeError:
  203.             self.instrument_name = None
  204.            
  205.         self.channels = []
  206.  
  207.         try:
  208.             conf = OrderedConfigParser()
  209.             conf.read(os.path.join(CONFIG_PATH, self.fullname+".cfg"))
  210.  
  211.             for section in conf.sections():
  212.                 if(not section[:-1].endswith("level") and
  213.                    not section.endswith("granules") and
  214.                    section.startswith(self.instrument_name)):
  215.                     name = eval(conf.get(section, "name"))
  216.                     try:
  217.                         w_range = eval(conf.get(section, "frequency"))
  218.                     except ConfigParser.NoOptionError:
  219.                         w_range = (-np.inf, -np.inf, -np.inf)
  220.                     try:
  221.                         resolution = eval(conf.get(section, "resolution"))
  222.                     except ConfigParser.NoOptionError:
  223.                         resolution = 0
  224.                     self.channels.append(Channel(name=name,
  225.                                                  wavelength_range=w_range,
  226.                                                  resolution=resolution))
  227.  
  228.         except (ConfigParser.NoSectionError, ConfigParser.NoOptionError):
  229.             for name, w_range, resolution in self.channel_list:
  230.                 self.channels.append(Channel(name=name,
  231.                                              wavelength_range=w_range,
  232.                                              resolution=resolution))
  233.  
  234.         self.channels_to_load = set([])
  235.  
  236.     def __getitem__(self, key, aslist = False):
  237.         if(isinstance(key, float)):
  238.             channels = [chn for chn in self.channels
  239.                         if(hasattr(chn, "wavelength_range") and
  240.                            chn.wavelength_range[0] <= key and
  241.                            chn.wavelength_range[2] >= key)]
  242.             channels = sorted(channels,
  243.                               lambda ch1,ch2:
  244.                                   ch1.__cmp__(ch2, key))
  245.            
  246.         elif(isinstance(key, str)):
  247.             channels = [chn for chn in self.channels
  248.                         if chn.name == key]
  249.             channels = sorted(channels)
  250.  
  251.         elif(isinstance(key, int)):
  252.             channels = [chn for chn in self.channels
  253.                         if int(np.round(chn.resolution)) == key]
  254.             channels = sorted(channels)
  255.  
  256.         elif(isinstance(key, (tuple, list))):
  257.             if len(key) == 0:
  258.                 raise KeyError("Key list must contain at least one element.")
  259.             channels = self.__getitem__(key[0], aslist = True)
  260.             if(len(key) > 1 and len(channels) > 0):
  261.                 dummy_instance = SatelliteInstrumentScene()
  262.                 dummy_instance.channels = channels
  263.                 channels = dummy_instance.__getitem__(key[1:], aslist = True)
  264.         else:
  265.             raise TypeError("Malformed key: " + str(key))
  266.  
  267.         if len(channels) == 0:
  268.             raise KeyError("No channel corresponding to "+str(key)+".")
  269.         elif aslist:
  270.             return channels
  271.         else:
  272.             return channels[0]
  273.  
  274.     def __setitem__(self, key, data):
  275.         # Add a channel if it is not already in the scene. Works only if key is
  276.         # a string.
  277.         try:
  278.             if key not in self:
  279.                 # if it's a blob with name and data, add it as is.
  280.                 if hasattr(data, "name") and hasattr(data, "data"):
  281.                     self.channels.append(data)
  282.                 else:
  283.                     kwargs = {"name": key}
  284.                     for attr in ["wavelength_range", "resolution"]:
  285.                         try:
  286.                             kwargs[attr] = getattr(data, attr)
  287.                         except (AttributeError, NameError):
  288.                             pass
  289.                     self.channels.append(Channel(**kwargs))
  290.         except AttributeError:
  291.             pass
  292.  
  293.         # Add the data.
  294.         if isinstance(data, np.ma.core.MaskedArray):
  295.             self[key].data = data
  296.         else:
  297.             try:
  298.                 self[key].data = data.data
  299.             except AttributeError:
  300.                 self[key].data = data
  301.  
  302.  
  303.  
  304.         # if isinstance(data, Channel):
  305.         #     self.channels.append(Channel(name=key,
  306.         #                              wavelength_range=data.wavelength_range,
  307.         #                              resolution=data.resolution))
  308.         #     self[key].data = data.data
  309.         # else:
  310.         #     try:
  311.         #         self[key].data = data
  312.         #     except KeyError:
  313.         #         self.channels.append(Channel(name=key))
  314.         #         self[key].data = data                            
  315.                
  316.  
  317.  
  318.     def __str__(self):
  319.         return "\n".join([str(chn) for chn in self.channels])
  320.  
  321.     def __iter__(self):
  322.         return self.channels.__iter__()
  323.  
  324.  
  325.     def _set_reader(self, pformat):
  326.         """Gets the reader for *pformat* format, and puts it in the `reader`
  327.        attribute.
  328.        """
  329.        
  330.         elements = pformat.split(".")
  331.         if len(elements) == 1:
  332.             reader_module = pformat
  333.             reader = "mpop.satin."+reader_module
  334.  
  335.             # Loading old style plugins
  336.             reader_module = pformat
  337.             LOG.info("old style plugin: " + pformat)
  338.             try:
  339.                 # Look for builtin reader
  340.                 loader = __import__(reader, globals(),
  341.                                     locals(), ['load'])
  342.             except ImportError:
  343.                 # Look for custom reader
  344.                 loader = __import__(reader_module, globals(),
  345.                                     locals(), ['load'])
  346.  
  347.             # Build a custom Reader plugin on the fly...
  348.             from mpop.plugin_base import Reader
  349.             reader_class = type(elements[-1].capitalize() + "Reader",
  350.                                 (Reader,),
  351.                                 {"pformat": elements[-1]})
  352.  
  353.             reader_instance = reader_class(self)
  354.  
  355.             # ... and set its "load" attribute with the "load" function of the
  356.             # loader module
  357.             loader = getattr(loader, "load")
  358.             setattr(reader_instance, "load", loader)
  359.  
  360.             setattr(self, elements[-1] + "_reader", reader_instance)
  361.  
  362.  
  363.         else:
  364.             reader_module = ".".join(elements[:-1])
  365.             reader_class = elements[-1]
  366.        
  367.             reader = "mpop.satin."+reader_module
  368.             try:
  369.                 # Look for builtin reader
  370.                 loader = __import__(reader, globals(),
  371.                                     locals(), [reader_class])
  372.             except ImportError:
  373.                 # Look for custom reader
  374.                 loader = __import__(reader_module, globals(),
  375.                                     locals(), [reader_class])
  376.             loader = getattr(loader, reader_class)
  377.             reader_instance = loader(self)
  378.             setattr(self, loader.pformat + "_reader", reader_instance)
  379.  
  380.         return reader_instance
  381.  
  382.  
  383.  
  384.  
  385.  
  386.     def load(self, channels=None, load_again=False, area_extent=None, **kwargs):
  387.         """Load instrument data into the *channels*. *Channels* is a list or a
  388.        tuple containing channels we will load data into, designated by there
  389.        center wavelength (float), resolution (integer) or name (string). If
  390.        None, all channels are loaded.
  391.  
  392.        The *load_again* boolean flag allows to reload the channels even they
  393.        have already been loaded, to mirror changes on disk for example. This
  394.        is false by default.
  395.  
  396.        The *area_extent* keyword lets you specify which part of the data to
  397.        load. Given as a 4-element sequence, it defines the area extent to load
  398.        in satellite projection.
  399.  
  400.        The other keyword arguments are passed as is to the reader
  401.        plugin. Check the corresponding documentation for more details.
  402.        """
  403.  
  404.         # Set up the list of channels to load.
  405.         if channels is None:
  406.             for chn in self.channels:
  407.                 self.channels_to_load |= set([chn.name])
  408.  
  409.         elif(isinstance(channels, (list, tuple, set))):
  410.             self.channels_to_load = set()
  411.             for chn in channels:
  412.                 try:
  413.                     self.channels_to_load |= set([self[chn].name])
  414.                 except KeyError:
  415.                     self.channels_to_load |= set([chn])
  416.  
  417.         else:
  418.             raise TypeError("Channels must be a list/"
  419.                             "tuple/set of channel keys!")
  420.  
  421.         loaded_channels = [chn.name for chn in self.loaded_channels()]
  422.         if load_again:
  423.             for chn in self.channels_to_load:
  424.                 if chn in loaded_channels:
  425.                     self.unload(chn)
  426.                     loaded_channels = []
  427.         else:
  428.             for chn in loaded_channels:
  429.                 self.channels_to_load -= set([chn])
  430.  
  431.         # find the plugin to use from the config file
  432.         conf = ConfigParser.ConfigParser()
  433.         try:
  434.             conf.read(os.path.join(CONFIG_PATH, self.fullname + ".cfg"))
  435.             if len(conf.sections()) == 0:
  436.                 raise ConfigParser.NoSectionError(("Config file did "
  437.                                                     "not make sense"))
  438.             levels = [section for section in conf.sections()
  439.                       if section.startswith(self.instrument_name+"-level")]
  440.         except ConfigParser.NoSectionError:
  441.             LOG.warning("Can't load data, no config file for " + self.fullname)
  442.             self.channels_to_load = set()
  443.             return
  444.        
  445.         levels.sort()
  446.  
  447.         if levels[0] == self.instrument_name+"-level1":
  448.             levels = levels[1:]
  449.  
  450.         if len(levels) == 0:
  451.             raise ConfigParser.NoSectionError(
  452.                 self.instrument_name+"-levelN (N>1) to tell me how to"+
  453.                 " read data... Not reading anything.")
  454.  
  455.         for level in levels:
  456.             if len(self.channels_to_load) == 0:
  457.                 return
  458.  
  459.             LOG.debug("Looking for sources in section "+level)
  460.             reader_name = conf.get(level, 'format')
  461.             try:
  462.                 reader_name = eval(reader_name)
  463.             except NameError:
  464.                 reader_name = str(reader_name)
  465.             LOG.debug("Using plugin mpop.satin."+reader_name)
  466.  
  467.             # read the data
  468.             reader = "mpop.satin."+reader_name
  469.             try:
  470.                 reader_instance = self._set_reader(reader_name)
  471.                 if area_extent is not None:
  472.                     if(isinstance(area_extent, (tuple, list)) and
  473.                        len(area_extent) == 4):
  474.                         kwargs["area_extent"] = area_extent
  475.                     else:
  476.                         raise ValueError("Area extent must be a sequence of "
  477.                                          "four numbers.")
  478.                    
  479.                 reader_instance.load(self, **kwargs)
  480.             except ImportError, e:
  481.                 LOG.exception("ImportError while loading "+reader_name+": "
  482.                               + str(e))
  483.                 continue
  484.             loaded_channels = set([chn.name for chn in self.loaded_channels()])
  485.             just_loaded = loaded_channels & self.channels_to_load
  486.             if len(just_loaded) == 0:
  487.                 LOG.info("No channels loaded with " + reader + ".")
  488.             self.channels_to_load -= loaded_channels
  489.             LOG.debug("Successfully loaded: "+str(just_loaded))
  490.            
  491.         if len(self.channels_to_load) > 0:
  492.             LOG.warning("Unable to import channels "
  493.                         + str(self.channels_to_load))
  494.  
  495.         self.channels_to_load = set()
  496.  
  497.     def save(self, filename, to_format="netcdf4", **options):
  498.         """Saves the current scene into a file of format *to_format*. Supported
  499.        formats are:
  500.        
  501.        - *netcdf4*: NetCDF4 with CF conventions.
  502.        """
  503.  
  504.         writer = "satout." + to_format
  505.         try:
  506.             writer_module = __import__(writer, globals(), locals(), ["save"])
  507.         except ImportError, err:
  508.             raise ImportError("Cannot load "+writer+" writer: "+str(err))
  509.  
  510.         return writer_module.save(self, filename, **options)
  511.  
  512.     def unload(self, *channels):
  513.         """Unloads *channels* from
  514.        memory. :meth:`mpop.scene.SatelliteInstrumentScene.load` must be called
  515.        again to reload the data.
  516.        """
  517.         for chn in channels:
  518.             try:
  519.                 self[chn].data = None
  520.             except AttributeError:
  521.                 LOG.warning("Can't unload channel" + str(chn))
  522.        
  523.        
  524.     def add_to_history(self, message):
  525.         """Adds a message to history info.
  526.        """
  527.         timestr = datetime.datetime.utcnow().isoformat()
  528.         timed_message = str(timestr + " - " + message)
  529.         if not self.info.get("history", ""):
  530.             self.info["history"] = timed_message
  531.         else:
  532.             self.info["history"] += "\n" + timed_message
  533.            
  534.  
  535.     def check_channels(self, *channels):
  536.         """Check if the *channels* are loaded, raise an error otherwise.
  537.        """
  538.         for chan in channels:
  539.             if not self[chan].is_loaded():
  540.                 raise NotLoadedError("Required channel %s not loaded,"
  541.                                      " aborting."%chan)
  542.  
  543.         return True
  544.  
  545.     def loaded_channels(self):
  546.         """Return the set of loaded_channels.
  547.        """
  548.         return set([chan for chan in self.channels if chan.is_loaded()])
  549.  
  550.     def project(self, dest_area, channels=None, precompute=False, mode=None, radius=None):
  551.         """Make a copy of the current snapshot projected onto the
  552.        *dest_area*. Available areas are defined in the region configuration
  553.        file (ACPG). *channels* tells which channels are to be projected, and
  554.        if None, all channels are projected and copied over to the return
  555.        snapshot.
  556.  
  557.        If *precompute* is set to true, the projecting data is saved on disk
  558.        for reusage. *mode* sets the mode to project in: 'quick' which works
  559.        between cartographic projections, and, as its denomination indicates,
  560.        is quick (but lower quality), and 'nearest' which uses nearest
  561.        neighbour for best projection. A *mode* set to None uses 'quick' when
  562.        possible, 'nearest' otherwise.
  563.  
  564.        *radius* defines the radius of influence for neighbour search in
  565.         'nearest' mode. Setting it to None, or omitting it will fallback to
  566.         default values (5 times the channel resolution) or 10km if the
  567.         resolution is not available.
  568.  
  569.        Note: channels have to be loaded to be projected, otherwise an
  570.        exception is raised.
  571.        """
  572.        
  573.         if not is_pyresample_loaded:
  574.             # Not much point in proceeding then
  575.             return self
  576.        
  577.         _channels = set([])
  578.  
  579.         if channels is None:
  580.             for chn in self.loaded_channels():
  581.                 _channels |= set([chn])
  582.  
  583.         elif(isinstance(channels, (list, tuple, set))):
  584.             for chn in channels:
  585.                 try:
  586.                     _channels |= set([self[chn]])
  587.                 except KeyError:
  588.                     LOG.warning("Channel "+str(chn)+" not found,"
  589.                                 "thus not projected.")
  590.         else:
  591.             raise TypeError("Channels must be a list/"
  592.                             "tuple/set of channel keys!")
  593.  
  594.  
  595.         res = copy.copy(self)
  596.  
  597.         if isinstance(dest_area, str):
  598.             dest_area = mpop.projector.get_area_def(dest_area)
  599.  
  600.        
  601.         res.area = dest_area
  602.         res.channels = []
  603.  
  604.         if not _channels <= self.loaded_channels():
  605.             LOG.warning("Cannot project nonloaded channels: %s."
  606.                         %(_channels - self.loaded_channels()))
  607.             LOG.info("Will project the other channels though.")
  608.             _channels = _channels and self.loaded_channels()
  609.        
  610.         cov = {}
  611.  
  612.         for chn in _channels:
  613.             if chn.area is None:
  614.                 if self.area is None:
  615.                     area_name = ("swath_" + self.fullname + "_" +
  616.                                  str(self.time_slot) + "_"
  617.                                  + str(chn.shape))
  618.                     chn.area = area_name
  619.                 else:
  620.                     if is_pyresample_loaded:
  621.                         try:                            
  622.                             chn.area = AreaDefinition(
  623.                                 self.area.area_id + str(chn.shape),
  624.                                 self.area.name,
  625.                                 self.area.proj_id,
  626.                                 self.area.proj_dict,
  627.                                 chn.shape[1],
  628.                                 chn.shape[0],
  629.                                 self.area.area_extent,
  630.                                 self.area.nprocs)
  631.    
  632.                         except AttributeError:
  633.                             try:
  634.                                 dummy = self.area.lons
  635.                                 dummy = self.area.lats
  636.                                 chn.area = self.area
  637.                                 area_name = ("swath_" + self.fullname + "_" +
  638.                                              str(self.time_slot) + "_"
  639.                                              + str(chn.shape))
  640.                                 chn.area.area_id = area_name
  641.                             except AttributeError:
  642.                                 chn.area = self.area + str(chn.shape)
  643.                     else:
  644.                         chn.area = self.area + str(chn.shape)
  645.             else: #chn.area is not None
  646.                 if (is_pyresample_loaded and
  647.                     isinstance(chn.area, SwathDefinition) and
  648.                     (not hasattr(chn.area, "area_id") or
  649.                      not chn.area.area_id)):
  650.                     area_name = ("swath_" + self.fullname + "_" +
  651.                                  str(self.time_slot) + "_"
  652.                                  + str(chn.shape) + "_"
  653.                                  + str(chn.name))
  654.                     chn.area.area_id = area_name
  655.  
  656.             if chn.area == dest_area:
  657.                 res.channels.append(chn)
  658.             else:
  659.                 if isinstance(chn.area, str):
  660.                     area_id = chn.area
  661.                 else:
  662.                     area_id = chn.area_id or chn.area.area_id
  663.                
  664.                 if area_id not in cov:
  665.                     if radius is None:
  666.                         if chn.resolution > 0:
  667.                             radius = 5 * chn.resolution
  668.                         else:
  669.                             radius = 10000
  670.                     cov[area_id] = mpop.projector.Projector(chn.area,
  671.                                                             dest_area,
  672.                                                             mode=mode,
  673.                                                             radius=radius)
  674.                     if precompute:
  675.                         try:
  676.                             cov[area_id].save()
  677.                         except IOError:
  678.                             LOG.exception("Could not save projection.")
  679.  
  680.                 try:
  681.                     res.channels.append(chn.project(cov[area_id]))
  682.                 except NotLoadedError:
  683.                     LOG.warning("Channel "+str(chn.name)+" not loaded, "
  684.                                 "thus not projected.")
  685.        
  686.         # Compose with image object
  687.         try:
  688.             if res._CompositerClass is not None:
  689.                 # Pass weak ref to compositor to allow garbage collection
  690.                 res.image = res._CompositerClass(weakref.proxy(res))
  691.         except AttributeError:
  692.             pass
  693.        
  694.         return res
  695.  
  696. if sys.version_info < (2, 5):
  697.     def any(iterable):
  698.         for element in iterable:
  699.             if element:
  700.                 return True
  701.         return False
  702.  
  703.  
  704. def assemble_segments(segments):
  705.     """Assemble the scene objects listed in *segment_list* and returns the
  706.    resulting scene object.
  707.    """
  708.     from mpop.satellites import GenericFactory
  709.    
  710.     channels = set([])
  711.     for seg in segments:
  712.         channels |= set([chn.name for chn in seg.loaded_channels()])
  713.  
  714.     seg = segments[0]
  715.    
  716.     new_scene = GenericFactory.create_scene(seg.satname, seg.number,
  717.                                             seg.instrument_name, seg.time_slot,
  718.                                             seg.orbit, variant=seg.variant)
  719.    
  720.     swath_definitions = {}
  721.  
  722.     for chn in channels:
  723.         new_scene[chn] = np.ma.concatenate([seg[chn].data
  724.                                             for seg in segments
  725.                                             if seg[chn].is_loaded()])
  726.         try:
  727.  
  728.             area_names = tuple([seg[chn].area.area_id
  729.                                 for seg in segments
  730.                                 if seg[chn].is_loaded()])
  731.             if area_names not in swath_definitions:
  732.            
  733.                 lons = np.ma.concatenate([seg[chn].area.lons[:]
  734.                                           for seg in segments
  735.                                           if seg[chn].is_loaded()])
  736.                 lats = np.ma.concatenate([seg[chn].area.lats[:]
  737.                                           for seg in segments
  738.                                           if seg[chn].is_loaded()])
  739.                 new_scene[chn].area = SwathDefinition(lons=lons, lats=lats)
  740.                 area_name = "+".join(area_names)
  741.                 new_scene[chn].area.area_id = area_name
  742.                 new_scene[chn].area_id = area_name
  743.                 swath_definitions[area_names] = new_scene[chn].area
  744.             else:
  745.                 new_scene[chn].area = swath_definitions[area_names]
  746.                 new_scene[chn].area_id = new_scene[chn].area.area_id
  747.         except AttributeError:
  748.             pass
  749.  
  750.     try:
  751.         lons = np.ma.concatenate([seg.area.lons[:] for seg in segments])
  752.         lats = np.ma.concatenate([seg.area.lats[:] for seg in segments])
  753.         new_scene.area = SwathDefinition(lons=lons, lats=lats)
  754.         for chn in channels:
  755.             if any([seg[chn].area for seg in segments]):
  756.                 try:
  757.                     lon_arrays = []
  758.                     lat_arrays = []
  759.                     for seg in segments:
  760.                         if seg[chn].area is not None:
  761.                             lon_arrays.append(seg[chn].area.lons[:])
  762.                             lat_arrays.append(seg[chn].area.lats[:])
  763.                         else:
  764.                             lon_arrays.append(seg.area.lons[:])
  765.                             lat_arrays.append(seg.area.lats[:])
  766.                     lons = np.ma.concatenate(lon_arrays)
  767.                     lats = np.ma.concatenate(lat_arrays)
  768.                     new_scene[chn].area = SwathDefinition(lons=lons, lats=lats)
  769.                 except AttributeError:
  770.                     pass
  771.     except AttributeError:
  772.         pass
  773.  
  774.     return new_scene
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement