Issues/35 (#38)

Removes warnings / errors when asked to interpolate the PWV within a gap of missing PWV data larger than a day.

pwv_kpno/__init__.py

@@ -53,18 +53,18 @@
 from . import package_settings
 from . import pwv_atm
 
-__authors__ = ['Daniel Perrefort']
-__copyright__ = 'Copyright 2017, Daniel Perrefort'
+__authors__ = ['MWV Research Group']
+__copyright__ = 'Copyright 2017, MWV Research Group'
 __credits__ = [
+    'Daniel Perrefort'
     'W. M. Wood-Vasey',
-    'K.Azalee Bostroem',
+    'K. Azalee Bostroem',
     'Jessica Kroboth',
     'Tom Matheson',
-    'Alexander Afanasyev',
 ]
 
 __license__ = 'GPL V3'
-__version__ = '1.2.0'
+__version__ = '1.2.1'
 __email__ = 'djperrefort@pitt.edu'
 __status__ = 'Release'
 

pwv_kpno/_download_pwv_data.py

@@ -30,14 +30,6 @@
 
 from .package_settings import settings
 
-__authors__ = ['Daniel Perrefort']
-__copyright__ = 'Copyright 2016, Daniel Perrefort'
-__credits__ = ['Jessica Kroboth']
-
-__license__ = 'GPL V3'
-__email__ = 'djperrefort@pitt.edu'
-__status__ = 'Release'
-
 
 @np.vectorize
 def _suomi_date_to_timestamp(year: int, days_str: str) -> float:

pwv_kpno/_update_pwv_model.py

@@ -33,15 +33,7 @@
 from ._download_pwv_data import update_local_data
 from .package_settings import settings
 
-__authors__ = ['Daniel Perrefort']
-__copyright__ = 'Copyright 2017, Daniel Perrefort'
-
-__license__ = 'GPL V3'
-__email__ = 'djperrefort@pitt.edu'
-__status__ = 'Release'
-
-warnings.filterwarnings("ignore",
-                        message='Empty data detected for ODR instance.')
+warnings.filterwarnings("ignore", message='Empty data detected for ODR instance.')
 
 
 def _linear_regression(x: np.array, y: np.array, sx: np.array, sy: np.array):

pwv_kpno/blackbody_with_atm.py

@@ -56,13 +56,6 @@
 
 from pwv_kpno.pwv_atm import trans_for_pwv
 
-__authors__ = ['Daniel Perrefort']
-__copyright__ = 'Copyright 2017, Daniel Perrefort'
-
-__license__ = 'GPL V3'
-__email__ = 'djperrefort@pitt.edu'
-__status__ = 'Release'
-
 
 def sed(temp: float,
         wavelengths: np.array,

pwv_kpno/exceptions.py

@@ -0,0 +1,24 @@
+#!/usr/bin/env python3
+# -*- coding: UTF-8 -*-
+
+#    This file is part of the pwv_kpno software package.
+#
+#    The pwv_kpno package 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.
+#
+#    The pwv_kpno package 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 pwv_kpno. If not, see <http://www.gnu.org/licenses/>.
+
+
+"""Custom exceptions for ``pwv_kpno``."""
+
+
+class ModelingConfigError(Exception):
+    pass

pwv_kpno/package_settings.py

@@ -72,19 +72,13 @@
 import os
 import shutil
 from datetime import datetime
-from warnings import warn, simplefilter
+from warnings import simplefilter, warn
 
 import numpy as np
+import scipy.interpolate as interpolate
 from astropy.table import Table
 
-from ._atm_model import create_pwv_atm_model
-
-__authors__ = ['Daniel Perrefort']
-__copyright__ = 'Copyright 2017, Daniel Perrefort'
-
-__license__ = 'GPL V3'
-__email__ = 'djperrefort@pitt.edu'
-__status__ = 'Release'
+from .exceptions import ModelingConfigError
 
 simplefilter('always', UserWarning)
 # Sites included with release that cannot be overwritten by the user
@@ -95,8 +89,60 @@
     'PWV', 'PWVerr', 'ZenithDelay', 'SrfcPress', 'SrfcTemp', 'SrfcRH')
 
 
-class ModelingConfigError(Exception):
-    pass
+def _calc_num_density_conversion():
+    """Calculate conversion factor from PWV * cross section to optical depth
+
+    Returns:
+        The conversion factor in units of 1 / (mm * cm^2)
+    """
+
+    n_a = 6.02214129E23  # 1 / mol (Avogadro's constant)
+    h2o_molar_mass = 18.0152  # g / mol
+    h2o_density = 0.99997  # g / cm^3
+    one_mm_in_cm = 10  # mm / cm
+
+    # Conversion factor 1 / (mm * cm^2)
+    mm_to_num_dens = (n_a * h2o_density) / (h2o_molar_mass * one_mm_in_cm)
+
+    return mm_to_num_dens
+
+
+def create_pwv_atm_model(
+        model_lambda: np.ndarray,
+        model_cs: np.ndarray,
+        out_lambda: np.ndarray) -> Table:
+    """Creates a table of conversion factors from PWV to optical depth
+
+    Expects input and output wavelengths to be in same units. Expects modeled
+    cross sections to be in cm^2.
+
+    Args:
+        model_lambda: Array of input wavelengths
+        model_cs: Array of cross sections for each input wavelength
+        out_lambda: Array of desired output wavelengths
+
+    Returns:
+        A table with columns 'wavelength' and '1/mm'
+    """
+
+    model_cs = np.array(model_cs)  # Just in case we are passed a list
+    if (model_cs < 0).any():
+        raise ValueError('Cross sections cannot be negative.')
+
+    if np.array_equal(model_lambda, out_lambda):
+        out_cs = model_cs  # This function requires ndarray behavior
+
+    else:
+        interp_cs = interpolate.interp1d(model_lambda, model_cs)
+        out_cs = interp_cs(out_lambda)
+
+    pwv_num_density = out_cs * _calc_num_density_conversion()
+    out_table = Table(
+        data=[out_lambda, pwv_num_density],
+        names=['wavelength', '1/mm']
+    )
+
+    return out_table
 
 
 def site_property(f):

pwv_kpno/pwv_atm.py

@@ -91,7 +91,7 @@
 """
 
 import os
-from datetime import datetime, timedelta
+from datetime import datetime
 from glob import glob
 from typing import List, Tuple, Union
 
@@ -105,55 +105,51 @@
 from ._update_pwv_model import update_models
 from .package_settings import settings
 
-__authors__ = ['Daniel Perrefort', 'Michael Wood-Vasey']
-__copyright__ = 'Copyright 2017, Daniel Perrefort'
-
-__license__ = 'GPL V3'
-__email__ = 'djperrefort@pitt.edu'
-__status__ = 'Release'
-
 
 def _warn_available_data(
-        test_dates: Union[float, np.array], known_dates: np.array):
-    """Check if a date falls within the range of data in an astropy table
+        test_dates: Union[float, np.array],
+        dates_with_data: np.array) -> None:
+    """Warn if given dates don't falls within a range dates with measurements
 
     Args:
-        test_dates: A timezone aware datetime
-        known_dates: An astropy table containing column 'date'
+        test_dates: Dates to check for available data
+        dates_with_data: Dates with data available
     """
 
     test_dates = np.atleast_1d(test_dates)  # In case passed a float
-    if not len(known_dates):
+    if not len(dates_with_data):
         err_msg = 'No PWV data for primary receiver available on local machine.'
         raise RuntimeError(err_msg)
 
     # Check date falls within the range of available PWV data
-    min_known_date, max_known_date = min(known_dates), max(known_dates)
-    dates_too_early = test_dates[test_dates < min_known_date]
-    if len(dates_too_early):
+    min_known_date = dates_with_data.min(),
+    if (test_dates < min_known_date).any():
         min_date = datetime.utcfromtimestamp(min_known_date)
         raise ValueError(
             f'No PWV data found for dates before {min_date} on local machine'
         )
 
-    dates_too_late = test_dates[test_dates > max_known_date]
-    if len(dates_too_late):
+    max_known_date = dates_with_data.max()
+    if (test_dates > max_known_date).any():
         max_date = datetime.utcfromtimestamp(max_known_date)
         raise ValueError(
             f'No PWV data found for dates after {max_date} on local machine'
         )
 
-    differences = (test_dates.reshape(1, -1) - known_dates.reshape(-1, 1))
-    indices = np.abs(differences).argmin(axis=0)
-    residual = np.diagonal(differences[indices,])
-
-    one_day_in_seconds = 24 * 60 * 60
-    out_of_interp_range = test_dates[residual > one_day_in_seconds]
-    if len(out_of_interp_range):
-        raise ValueError(
-            f'Specified datetimes falls within interval of missing SuomiNet' 
-            f' data larger than 1 day: {out_of_interp_range}.'
-        )
+    # Todo: Warn if dat falls in interval > 1 day
+    # The code below is extremely slow
+    # one_day_in_seconds = 24 * 60 * 60
+    # interval_in_seconds = interval * one_day_in_seconds
+    # out_of_interp_range = []
+    # for date in test_dates:
+    #     if (np.abs(dates_with_data - date) > interval_in_seconds).any():
+    #         out_of_interp_range.append(date)
+    #
+    # if out_of_interp_range:
+    #     raise ValueError(
+    #         f'Specified datetimes falls within interval of missing SuomiNet'
+    #         f' data larger than 1 day: {out_of_interp_range}.'
+    #     )
 
 
 def _pwv_date(
@@ -522,7 +518,7 @@ def _trans_for_date(
 
 def trans_for_date(
         date: Union[float, np.array, datetime],
-        airmass: float = 1, 
+        airmass: float = 1,
         format: str = None,
         bins: Union[int, list] = None) -> Table:
     """Return a model for the atmospheric transmission function due to PWV

tests/_create_mock_data.py

@@ -23,13 +23,6 @@
 from astropy.table import Table
 from pytz import utc
 
-__authors__ = ['Daniel Perrefort']
-__copyright__ = 'Copyright 2017, Daniel Perrefort'
-
-__license__ = 'GPL V3'
-__email__ = 'djperrefort@pitt.edu'
-__status__ = 'Release'
-
 
 def create_mock_pwv_model(year, gaps=None):
     """Create a mock model for the PWV level at Kitt Peak for airmass 1