added get_time function to calculate time for a given solar position

Author: alorenzo175

.gitignore

@@ -4,6 +4,7 @@ __pycache__/
 
 # emacs temp files
 *~
+#*#
 
 # C extensions
 *.so
@@ -68,6 +69,6 @@ coverage.xml
 # vi
 *.swp
 
-#Ignore notebooks
+#Ignore some notebooks
 *.ipynb
 !docs/tutorials/*.ipynb

LICENSE

@@ -74,7 +74,7 @@ import pvlib.solarposition
 import pvlib.clearsky
 
 # make a location
-tus = Location(32.2, -111, 700, 'MST')
+tus = Location(32.2, -111, 'MST', 700)
 
 # make a pandas DatetimeIndex for some day
 times = pd.date_range(start=datetime.datetime(2014,6,24), end=datetime.datetime(2014,6,25), freq='1Min')

MANIFEST.in

@@ -4,11 +4,13 @@
 import logging
 pvl_logger = logging.getLogger('pvlib')
 
+import datetime as dt
 import pdb
 import ast
 import re
 
 import numpy as np 
+import pytz
 
 
 class repack():   #repack a dict as a struct
@@ -303,4 +305,86 @@ def asind(number):
     return res
 
 
+def localize_to_utc(time, location):
+    """
+    Converts or localizes a time series to UTC.
+    
+    Parameters
+    ----------
+    time : datetime.datetime, pandas.DatetimeIndex, 
+           or pandas.Series/DataFrame with a DatetimeIndex.
+    location : pvlib.Location object
+    
+    Returns
+    -------
+    pandas object localized to UTC.
+    """
+    import datetime as dt
+    import pytz
+
+    if isinstance(time, dt.datetime):
+        if time.tzinfo is None:
+            time = pytz.timezone(location.tz).localize(time)
+        time_utc = time.astimezone(pytz.utc)
+    else:
+        try:
+            time_utc = time.tz_convert('UTC')
+            pvl_logger.debug('tz_convert to UTC')
+        except TypeError:
+            time_utc = time.tz_localize(location.tz).tz_convert('UTC')
+            pvl_logger.debug('tz_localize to {} and then tz_convert to UTC'
+                             .format(location.tz))
+        
+        
+    return time_utc
+
+
+def datetime_to_djd(time):
+    """
+    Converts a datetime to the Dublin Julian Day
+
+    Parameters
+    ----------
+    time : datetime.datetime 
+        time to convert
+
+    Returns
+    -------
+    float 
+        fractional days since 12/31/1899+0000
+    """
+
+    if time.tzinfo is None:
+        time_utc = pytz.utc.localize(time)
+    else:
+        time_utc = time.astimezone(pytz.utc)
+        
+    djd_start = pytz.utc.localize(dt.datetime(1899, 12, 31, 12))
+    djd = (time_utc - djd_start).total_seconds() * 1.0/(60 * 60 * 24)
+
+    return djd
+
+
+def djd_to_datetime(djd, tz='UTC'):
+    """
+    Converts a Dublin Julian Day float to a datetime.datetime object
+
+    Parameters
+    ----------
+    djd : float
+        fractional days since 12/31/1899+0000
+    tz : str
+        timezone to localize the result to
+
+    Returns
+    -------
+    datetime.datetime 
+       The resultant datetime localized to tz
+    """
+    
+    djd_start = pytz.utc.localize(dt.datetime(1899, 12, 31, 12))
 
+    utc_time = djd_start + dt.timedelta(days=djd)
+    return utc_time.astimezone(pytz.timezone(tz))
+    
+    

README.md

@@ -7,29 +7,20 @@
 # Will Holmgren (@wholmgren), University of Arizona, 2014 
 
 from __future__ import division
-
 import logging
 pvl_logger = logging.getLogger('pvlib')
+import datetime as dt
 
-import datetime
 
 import numpy as np
 import pandas as pd
 import pytz
 
-try:
-    from .spa_c_files.spa_py import spa_calc
-except ImportError as e:
-    pvl_logger.exception('Could not import built-in SPA calculator. You may need to recompile the SPA code.')
-
-try:    
-    import ephem
-except ImportError as e:
-    pvl_logger.warning('PyEphem not found.')
 
+from .pvl_tools import localize_to_utc, datetime_to_djd, djd_to_datetime
 
 
-def get_solarposition(time, location, method='spa', pressure=101325, 
+def get_solarposition(time, location, method='pyephem', pressure=101325, 
                       temperature=12):
     """
     A convenience wrapper for the solar position calculators.
@@ -48,6 +39,8 @@ def get_solarposition(time, location, method='spa', pressure=101325,
     """
 
     method = method.lower()
+    if isinstance(time, dt.datetime):
+        time = pd.DatetimeIndex([time,])
 
     if method == 'spa':
         ephem_df = spa(time, location)
@@ -59,7 +52,6 @@ def get_solarposition(time, location, method='spa', pressure=101325,
     return ephem_df
 
 
-
 def spa(time, location, raw_spa_output=False):
     '''
     Calculate the solar position using the C implementation of the NREL 
@@ -90,10 +82,16 @@ def spa(time, location, raw_spa_output=False):
 
     # Added by Rob Andrews (@Calama-Consulting), Calama Consulting, 2014 
     # Edited by Will Holmgren (@wholmgren), University of Arizona, 2014 
+
+    try:
+        from .spa_c_files.spa_py import spa_calc
+    except ImportError as e:
+        raise ImportError('Could not import built-in SPA calculator. '+
+                          'You may need to recompile the SPA code.')
 
     pvl_logger.debug('using built-in spa code to calculate solar position')
 
-    time_utc = _localize_to_utc(time, location)
+    time_utc = localize_to_utc(time, location)
 
     spa_out = []
 
@@ -120,6 +118,20 @@ def spa(time, location, raw_spa_output=False):
         return dfout
 
 
+def _ephem_setup(location, pressure, temperature):
+    import ephem
+    # initialize a PyEphem observer
+    obs = ephem.Observer()
+    obs.lat = str(location.latitude)
+    obs.lon = str(location.longitude)
+    obs.elevation = location.altitude
+    obs.pressure = pressure / 100. # convert to mBar
+    obs.temp = temperature
+    
+    # the PyEphem sun
+    sun = ephem.Sun()
+    return obs, sun
+
 
 def pyephem(time, location, pressure=101325, temperature=12):
     """
@@ -144,24 +156,17 @@ def pyephem(time, location, pressure=101325, temperature=12):
     """
 
     # Written by Will Holmgren (@wholmgren), University of Arizona, 2014 
-    
+
+    import ephem
+
     pvl_logger.debug('using PyEphem to calculate solar position')
 
-    time_utc = _localize_to_utc(time, location)
+    time_utc = localize_to_utc(time, location)
 
     sun_coords = pd.DataFrame(index=time_utc)
 
-    # initialize a PyEphem observer
-    obs = ephem.Observer()
-    obs.lat = str(location.latitude)
-    obs.lon = str(location.longitude)
-    obs.elevation = location.altitude
-    obs.pressure = pressure / 100. # convert to mBar
-    obs.temp = temperature
-    
-    # the PyEphem sun
-    sun = ephem.Sun()
-    
+    obs, sun = _ephem_setup(location, pressure, temperature)
+
     # make and fill lists of the sun's altitude and azimuth
     # this is the pressure and temperature corrected apparent alt/az.
     alts = []
@@ -197,8 +202,7 @@ def pyephem(time, location, pressure=101325, temperature=12):
         return sun_coords.tz_convert(location.tz)
     except TypeError:
         return sun_coords.tz_localize(location.tz)
-        
-        
+
 
 def ephemeris(time, location, pressure=101325, temperature=12):
     ''' 
@@ -217,7 +221,7 @@ def ephemeris(time, location, pressure=101325, temperature=12):
     pressure : float or DataFrame
           Ambient pressure (Pascals)
 
-    tempreature: float or DataFrame
+    temperature : float or DataFrame
           Ambient temperature (C)
       
     Returns
@@ -374,30 +378,63 @@ def ephemeris(time, location, pressure=101325, temperature=12):
     DFOut['solar_time'] = SolarTime
 
     return DFOut
-        
-        
-        
-def _localize_to_utc(time, location):
+
+
+def calc_time(lower_bound, upper_bound, location, attribute, value, 
+              pressure=101325, temperature=12, xtol=1.0e-12):
     """
-    Converts or localizes a time series to UTC.
+    Calculate the time between lower_bound and upper_bound
+    where the attribute is equal to value. Uses PyEphem for
+    solar position calculations.
     
     Parameters
     ----------
-    time : pandas.DatetimeIndex or pandas.Series/DataFrame with a DatetimeIndex.
+    lower_bound : datetime.datetime 
+    upper_bound : datetime.datetime 
     location : pvlib.Location object
+    attribute : str
+        The attribute of a pyephem.Sun object that
+        you want to solve for. Likely options are 'alt'
+        and 'az' (which must be given in radians).
+    value : int or float
+        The value of the attribute to solve for
+    pressure : int or float, optional
+        Air pressure in Pascals. Set to 0 for no
+        atmospheric correction.
+    temperature : int or float, optional
+        Air temperature in degrees C.
+    xtol : float, optional
+        The allowed error in the result from value
     
     Returns
     -------
-    pandas object localized to UTC.
+    datetime.datetime
+
+    Raises
+    ------
+    ValueError
+        If the value is not contained between the bounds.
+    AttributeError
+        If the given attribute is not an attribute of a
+        PyEphem.Sun object.
     """
-    
+
     try:
-        time_utc = time.tz_convert('UTC')
-        pvl_logger.debug('tz_convert to UTC')
-    except TypeError:
-        time_utc = time.tz_localize(location.tz).tz_convert('UTC')
-        pvl_logger.debug('tz_localize to {} and then tz_convert to UTC'
-                         .format(location.tz))
+        import scipy.optimize as so
+    except ImportError as e:
+        raise ImportError('The calc_time function requires scipy')
+
+    obs, sun = _ephem_setup(location, pressure, temperature)
+
+    def compute_attr(thetime, target, attr):
+        obs.date = thetime
+        sun.compute(obs)
+        return getattr(sun, attr) - target
 
-    return time_utc
-    
+    lb = datetime_to_djd(lower_bound)
+    ub = datetime_to_djd(upper_bound)
+
+    djd_root = so.brentq(compute_attr, lb, ub, 
+                         (value, attribute), xtol=xtol)
+
+    return djd_to_datetime(djd_root, location.tz)

pvlib/__init__.py

@@ -76,6 +76,28 @@ def test_pyephem_physical_dst():
 def test_pyephem_localization():  
     assert_frame_equal(solarposition.pyephem(times, tus), solarposition.pyephem(times_localized, tus))
 
+
+def test_calc_time():
+    import pytz
+    import math
+    # validation from USNO solar position calculator online
+
+    epoch = datetime.datetime(1970,1,1)
+    epoch_dt = pytz.utc.localize(epoch)
+    
+    loc = tus
+    loc.pressure = 0
+    actual_time = pytz.timezone(loc.tz).localize(datetime.datetime(2014, 10, 10, 8, 30))
+    lb = pytz.timezone(loc.tz).localize(datetime.datetime(2014, 10, 10, 6))
+    ub = pytz.timezone(loc.tz).localize(datetime.datetime(2014, 10, 10, 10))
+    alt = solarposition.calc_time(lb, ub, loc, 'alt', math.radians(24.7))
+    az = solarposition.calc_time(lb, ub, loc, 'az', math.radians(116.3))
+    actual_timestamp = (actual_time - epoch_dt).total_seconds()
+    
+    assert_almost_equals((alt.replace(second=0, microsecond=0) - 
+                          epoch_dt).total_seconds(), actual_timestamp)
+    assert_almost_equals((az.replace(second=0, microsecond=0) - 
+                          epoch_dt).total_seconds(), actual_timestamp)
+    
 
-        
-# add tests for daylight savings time?
+# add tests for daylight savings time?