Merge branch 'pre-master' of github.com:mraspaud/mpop into pre-master

mpop/satin/aapp1b.py

@@ -10,6 +10,7 @@
 #   Nina Håkansson <nina.hakansson@smhi.se>
 #   Oana Nicola <oananicola@yahoo.com>
 #   Lars Ørum Rasmussen <ras@dmi.dk>
+#   Panu Lahtinen <panu.lahtinen@fmi.fi>
 
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -39,7 +40,11 @@
 from ConfigParser import ConfigParser
 from mpop import CONFIG_PATH
 
+<<<<<<< HEAD
+LOGGER = logging.getLogger('aapp1b')
+=======
 logger = logging.getLogger('aapp1b')
+>>>>>>> upstream/pre-master
 
 def load(satscene, *args, **kwargs):
     """Read data from file and load it into *satscene*.
@@ -58,7 +63,11 @@ def load(satscene, *args, **kwargs):
 
     options["calibrate"] = kwargs.get("calibrate", True)
 
+<<<<<<< HEAD
+    LOGGER.info("Loading instrument '%s'" % satscene.instrument_name)
+=======
     logger.info("Loading instrument '%s'" % satscene.instrument_name)
+>>>>>>> upstream/pre-master
     try:
         CASES[satscene.instrument_name](satscene, options)
     except KeyError:
@@ -71,7 +80,11 @@ def load_avhrr(satscene, options):
         raise IOError("No filename given, cannot load.")
 
     chns = satscene.channels_to_load & set(AVHRR_CHANNEL_NAMES)
+<<<<<<< HEAD
+    LOGGER.info("Loading channels " + str(sorted(list(chns))))
+=======
     logger.info("Loading channels " + str(sorted(list(chns))))
+>>>>>>> upstream/pre-master
     if len(chns) == 0:
         return
 
@@ -96,7 +109,11 @@ def load_avhrr(satscene, options):
 
     filename = file_list[0]
 
+<<<<<<< HEAD
+    LOGGER.debug("Loading from " + filename)
+=======
     logger.debug("Loading from " + filename)
+>>>>>>> upstream/pre-master
 
     scene = AAPP1b(filename)
     scene.read()
@@ -107,7 +124,11 @@ def load_avhrr(satscene, options):
     try:
         from pyresample import geometry
     except ImportError, ex_:
+<<<<<<< HEAD
+        LOGGER.debug("Could not load pyresample: " + str(ex_))
+=======
         logger.debug("Could not load pyresample: " + str(ex_))
+>>>>>>> upstream/pre-master
         satscene.lat = scene.lats
         satscene.lon = scene.lons
     else:
@@ -119,6 +140,7 @@ def load_avhrr(satscene, options):
             satscene[chn].info['units'] = scene.units[chn]
             satscene[chn].area = satscene.area
 
+
 AVHRR_CHANNEL_NAMES = ("1", "2", "3A", "3B", "4", "5")
 
 # AAPP 1b header
@@ -321,7 +343,11 @@ def read(self):
             fp_.seek(10664 * 2, 1)
             data = np.fromfile(fp_, dtype=_SCANTYPE)
 
+<<<<<<< HEAD
+        LOGGER.debug("Reading time " + str(datetime.datetime.now() - tic))
+=======
         logger.debug("Reading time " + str(datetime.datetime.now() - tic))
+>>>>>>> upstream/pre-master
         self._header = header
         self._data = data
 
@@ -335,7 +361,11 @@ def navigate(self):
         try:
             from geotiepoints import SatelliteInterpolator
         except ImportError:
+<<<<<<< HEAD
+            LOGGER.warning("Could not interpolate lon/lats, "
+=======
             logger.warning("Could not interpolate lon/lats, "
+>>>>>>> upstream/pre-master
                         "python-geotiepoints missing.")
             self.lons, self.lats = lons40km, lats40km
         else:
@@ -354,7 +384,11 @@ def navigate(self):
                                            along_track_order,
                                            cross_track_order)
             self.lons, self.lats = satint.interpolate()
+<<<<<<< HEAD
+            LOGGER.debug("Navigation time " + str(datetime.datetime.now() - tic))
+=======
             logger.debug("Navigation time " + str(datetime.datetime.now() - tic))
+>>>>>>> upstream/pre-master
 
     def calibrate(self, chns=("1", "2", "3A", "3B", "4", "5"),
                   calibrate=1):
@@ -424,7 +458,11 @@ def calibrate(self, chns=("1", "2", "3A", "3B", "4", "5"),
             else:
                 self.units['5'] = ''
 
+<<<<<<< HEAD
+        LOGGER.debug("Calibration time " + str(datetime.datetime.now() - tic))
+=======
         logger.debug("Calibration time " + str(datetime.datetime.now() - tic))
+>>>>>>> upstream/pre-master
 
 
 def _vis_calibrate(data, chn, calib_type):
@@ -441,7 +479,11 @@ def _vis_calibrate(data, chn, calib_type):
         return channel
 
     if calib_type == 2:
+<<<<<<< HEAD
+        LOGGER.info("Radiances are not yet supported for " + 
+=======
         logger.info("Radiances are not yet supported for " + 
+>>>>>>> upstream/pre-master
                  "the VIS/NIR channels!")
 
     mask1 = channel <= np.expand_dims(data["calvis"][:, chn, 2, 4], 1)
@@ -503,8 +545,24 @@ def _ir_calibrate(header, data, irchn, calib_type):
     bandcor_2 = header['radtempcnv'][0, irchn, 1]/1e5
     bandcor_3 = header['radtempcnv'][0, irchn, 2]/1e6
 
+<<<<<<< HEAD
+    # Count to radiance conversion:
+    rad = k1_ * count*count + k2_*count + k3_
+
+    if calib_type == 2:
+        return rad
+
+    all_zero = np.logical_and(np.logical_and(np.equal(k1_, 0),
+                                             np.equal(k2_, 0)),
+                              np.equal(k3_, 0))    
+    idx = np.indices((all_zero.shape[0],))
+    suspect_line_nums = np.repeat(idx[0], all_zero[:, 0])
+    if suspect_line_nums.any():
+        LOGGER.info("Suspect scan lines: " + str(suspect_line_nums))
+=======
     ir_const_1 = 1.1910659e-5
     ir_const_2 = 1.438833
+>>>>>>> upstream/pre-master
 
     t_planck = (ir_const_2*cwnum) / np.log(1 + ir_const_1*cwnum*cwnum*cwnum/rad)
 

mpop/satin/helper_functions.py

@@ -0,0 +1,126 @@
+
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2014.
+# 
+# Author(s):
+#  
+#   Panu Lahtinen <panu.lahtinen@fmi.fi>
+#
+# 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/>.
+
+'''Helper functions for area extent calculations.
+'''
+
+import numpy as np
+from mpop.projector import get_area_def
+# from pyresample.utils import AreaNotFound
+
+import pyresample
+
+import logging
+from pyproj import Proj
+
+LOGGER = logging.getLogger(__name__)
+
+
+def area_def_names_to_extent(area_def_names, proj4_str,
+                             default_extent=(-5567248.07, -5570248.48,
+                                              5570248.48, 5567248.07)):
+    '''Convert a list of *area_def_names* to maximal area extent in
+    destination projection defined by *proj4_str*. *default_extent*
+    gives the extreme values.  Default value is MSG3 extents at
+    lat0=0.0.
+    '''
+
+    if type(area_def_names) is not list:
+        area_def_names = [area_def_names]
+
+    maximum_extent = None
+
+    for name in area_def_names:
+
+        try:
+            boundaries = get_area_def(name).get_boundary_lonlats()
+        except pyresample.utils.AreaNotFound:
+            LOGGER.warning('Area definition not found ' + name)
+            continue
+        except AttributeError:
+            boundaries = name.get_boundary_lonlats()
+
+        lon_sides = (boundaries[0].side1, boundaries[0].side2,
+                     boundaries[0].side3, boundaries[0].side4)
+        lat_sides = (boundaries[1].side1, boundaries[1].side2,
+                     boundaries[1].side3, boundaries[1].side4)
+
+
+        maximum_extent = boundaries_to_extent(proj4_str, maximum_extent,
+                                              default_extent,
+                                              lon_sides, lat_sides)
+
+    return maximum_extent
+
+
+def boundaries_to_extent(proj4_str, maximum_extent, default_extent,
+                         lon_sides, lat_sides):
+    '''Get area extent from given boundaries.
+    '''
+
+    # proj4-ify the projection string
+    if '+' not in proj4_str:
+        proj4_str = proj4_str.split(' ')
+        proj4_str = '+' + ' +'.join(proj4_str)
+
+    pro = Proj(proj4_str)
+
+    # extents for edges
+    x_dir, y_dir = pro(np.concatenate(lon_sides),
+                       np.concatenate(lat_sides))
+
+    # replace invalid values with NaN
+    x_dir[np.abs(x_dir) > 1e20] = np.nan
+    y_dir[np.abs(y_dir) > 1e20] = np.nan
+
+    # Get the maximum needed extent from different corners.
+    extent = [np.nanmin(x_dir),
+              np.nanmin(y_dir),
+              np.nanmax(x_dir),
+              np.nanmax(y_dir)]
+
+    # Replace "infinity" values with default extent
+    for i in range(4):
+        if extent[i] is np.nan:
+            extent[i] = default_extent[i]
+
+    # update maximum extent
+    if maximum_extent is None:
+        maximum_extent = extent
+    else:
+        if maximum_extent[0] > extent[0]:
+            maximum_extent[0] = extent[0]
+        if maximum_extent[1] > extent[1]:
+            maximum_extent[1] = extent[1]
+        if maximum_extent[2] < extent[2]:
+            maximum_extent[2] = extent[2]
+        if maximum_extent[3] < extent[3]:
+            maximum_extent[3] = extent[3]
+
+    # Replace "infinity" values with default extent
+    for i in range(4):
+        if not np.isfinite(maximum_extent[i]):
+            maximum_extent[i] = default_extent[i]
+
+    return maximum_extent
+