Skip to content
Snippets Groups Projects
Commit 05d38736 authored by Bruce Flynn's avatar Bruce Flynn
Browse files

No commit message

No commit message
parents
Branches
Tags
No related merge requests found
Makefile 0 → 100644
ifndef PYBIN
PYBIN=python
endif
all: clean egg
egg:
$(PYBIN) setup.py bdist_egg
clean:
rm -rf dist
rm -rf build
#rm -rf *.egg-info
torepos:
scp dist/*.egg $$USER@bora7.ssec.wisc.edu:/home/httpd/html/eggs/repos/MetObs
develop:
$(PYBIN) setup.py develop
undevelop:
$(PYBIN) setup.py develop --uninstall
__import__('pkg_resources').declare_namespace(__name__)
\ No newline at end of file
"""Functions for calculating the dervied data.
"""
__version__ = '$Revision: 1.15 $'
#$Source: /cvsroot/TOOLS/dev/metobs/python/data/metobs/data/__init__.py,v $
__docformat__ = 'restructuredtext en'
import math
from datetime import timedelta
import fpconst
from numpy import array, average, zeros, radians, degrees, arctan2, sin, cos
from pydap.client import open_url as dapopen
from metobs import mytime
#: Not a number representation
NaN = fpconst.NaN
#: function to determine if a number is NaN
is_nan = fpconst.isNaN
def average_for_interval(basetime, arr, interval):
"""Generate averages for timestamped data at a particular interval. Averages
may be produced from different numbers of points depending on the number of
points in each interval bucket. Returns data averaged at the specified
interval. Times are averaged as well as the data meaning the time will
represent the middle of the average.
:param basetime: Time to use as the base for computing intervals
:type basetime: ``datetime.datetime``
:param arr: array of data where index 0 in each row is an epoch time in seconds.
:type arr: ``numpy.array``
:param interval: the averaging interval in seconds
"""
arr_out = []
stop_dt = basetime + timedelta(seconds=interval)
stop_idx = 1
start_idx = 0
tz_cls = stop_dt.tzinfo and stop_dt.tzinfo.__class__
# advance through data until there is a timestamp in the averaging interval
while mytime.seconds_to_datetime(arr[0][0], tz=tz_cls) >= stop_dt:
stop_dt += timedelta(seconds=interval)
for row_idx in range(arr.shape[0]):
row = arr[row_idx]
cur_dt = mytime.seconds_to_datetime(row[0], tz=tz_cls)
# determine stop index
if cur_dt < stop_dt:
stop_idx +=1
# stop and average
else:
s = mytime.seconds_to_datetime(arr[start_idx][0], tz=tz_cls)
# matrix of data to be averaged
avg_input = arr[start_idx:stop_idx,1:]
avg_out = array(zeros(arr[0].shape[0]))
avg_out.fill(NaN)
# timestamp for center of averaging interval
int_dt = stop_dt - timedelta(seconds=(interval/2))
avg_out[0] = mytime.datetime_to_epoch(int_dt)
avg_out[1:] = average(avg_input, axis=0)
arr_out.append(avg_out)
# adjust conditions for next average
start_idx = row_idx
stop_idx = start_idx
stop_dt += timedelta(seconds=interval)
avg_arr = arr[start_idx:,:]
arr_out.append(average(avg_arr, axis=0))
return array(arr_out)
def average_for_interval_degree(basetime, arr, interval):
"""Generate averages for timestamped data at a particular interval.
Data must be in degrees. See documentation for average_for_interval()
for more information on the averaging of this function.
:param basetime: Time to use as the base for computing intervals
:type basetime: ``datetime.datetime``
:param arr: array of data in degrees where index 0 in each row is an epoch time in seconds.
:type arr: ``numpy.array``
:param interval: the averaging interval in seconds
"""
# Create 2 x and y component arrays from radians
y_comp = array((arr[:,0],sin(radians(arr[:,1])))).transpose()
x_comp = array((arr[:,0],cos(radians(arr[:,1])))).transpose()
# Average the 2 components
y_comp = average_for_interval(basetime, y_comp, interval)
x_comp = average_for_interval(basetime, x_comp, interval)
# Convert back to degrees and then add 360 to negative degrees
new_degree = array((y_comp[:,0],degrees(arctan2(y_comp[:,1], x_comp[:,1]))))
for i in range(new_degree[1].size):
if new_degree[1,i] < 0:
new_degree[1,i] += 360
new_degree = new_degree.transpose()
return new_degree
def get_nc_array(url, vars):
"""Get an array of data for a date and set of vars.
:param url: where to retreive data from
:param vars: netcdf variable names
:type vars: ``list``
:raise ValueError: If variable name is not in file
:raise IndexError: If variable has no data
"""
nc = dapopen(url)
arr = []
for nm in vars:
if nm not in nc.keys():
raise ValueError("'%s' not a variable in %s" %(nm, nc.name))
if nc[nm].shape[0] == 0:
raise IndexError("No data for '%s'" % nm)
arr.append(nc[nm][:])
return array(arr).transpose()
def approx_eq(this, that, delta=.000001):
"""Check if to decimal numbers are approximately equal.
:param this: operand one
:param that: operand two
:param delta: amount of allowable difference
"""
if this > that - delta and this < that + delta:
return True
return False
def dewpoint(tempC, relhum):
"""
Algorithm from Tom Whittaker tempC is the temperature in degrees Celsius,
relhum is the relative humidity as a percentage.
:param tempC: temperature in celsius
:param relhum: relative humidity as a percentage
"""
if not tempC or not relhum:
return NaN
gasconst = 461.5
latheat = 2500800.0
dp = 1.0 / ( 1.0 / ( 273.15 + tempC ) - gasconst * math.log( (0.0 + relhum) / 100 ) / \
( latheat - tempC * 2397.5 ))
return min(dp - 273.15, tempC)
def relhum(airTempK, dewpointTempK):
"""
Algorithm from derived by David Hoese from the above
dewpoint(tempC, relhum) function, both parameters are in Kelvin units.
:param airTempK: air temperature in Kelvin
:param dewpointTempK: dewpoint temp in Kelvin
"""
if airTempK == None or dewpointTempK == None:
return NaN
gas_constant = 461.5
latheat = 2500800.0
# Only one section of the equation
latpart = (latheat-(airTempK-273.15)*2397.5)
relativehum = 100 * math.e**((latpart/airTempK - latpart/dewpointTempK)/gas_constant)
return relativehum
def potentialtemp(airTempK, pressureMB):
"""
Algorithm from David Hoese to calculate potential temperature.
:param airTempK: air temperature in Kelvin
:param pressureMB: air pressure in millibars
"""
if airTempK == None or pressureMB == None:
return NaN
pT = airTempK * (pressure.max()/pressure)**.286
return pT
def awips_dewpoint(tempC, relhum):
"""
Algorithm taken from http://meted.ucar.edu/awips/validate/dewpnt.htm
tempC is the temperature in degrees celsius, relhum is the relative humidity
as a percentage.
:param tempC: temperature in celsius
:param relhum: relative humidity as a percentage
"""
if not tempC or not relhum:
return NaN
C15 = 26.66082
C1 = 0.0091379024
C2 = 6106.396
C3 = 223.1986
C4 = 0.0182758048
t = tempC + 273.15 # convert temperature to Kelvin
rh = relhum / 100 # convert relative humidity to ratio
es = math.exp( C15 - C1 * t - C2 / t ) # saturation vapor pressure
e = rh * es
b = C15 - math.log(e)
td = ( b - math.sqrt( b * b - C3 ) ) / C4
return min(td - 273.15, tempC)
def linearConvert(data, factor=1, offset=0):
"""
>>> linearConvert(0)
0.0
>>> linearConvert(1)
1.0
>>> linearConvert(1, 1, 0)
1.0
>>> linearConvert(1, 90, 10)
100.0
"""
return float(data) * factor + offset
def c2f (val):
"""Degrees celsius to fahrenheit.
>>> c2f(0)
32.0
>>> c2f(100)
212.0
"""
return linearConvert(val, (9/5.), 32)
def mps2mph (val):
"""Speed in meters per second to miles per hour.
"""
return linearConvert(val, 2.23694)
def mps2knots (val):
"""Meters per second to knots.
"""
return linearConvert(val, 1.9438445)
def altimeter(val, alt):
"""Compute altimeter from pressure and altitude.
:param val: pressure in hPa.
:param alt: altitude of the instrument in meters.
:returns: altimeter in inHg
"""
n = .190284
H = 327.538
c1= .0065 * pow(1013.25, n)/288.
c2 = H / pow( (p-.3), n)
ff = pow( 1.+c1*c2, 1./n)
return ((p-.3) * ff * 29.92 / 1013.25)
def dir2txt (val):
"""Convert degrees [0, 360) to a textual representation.
:param val: decimal degrees
>>> dir2txt(0)
'N'
>>> dir2txt(90)
'E'
>>> dir2txt(180)
'S'
>>> dir2txt(270)
'W'
>>> dir2txt(359)
'N'
"""
assert val >= 0 and val < 360, "'%s' out of range" % val
dirs = ("NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW")
if ((val >= 348.75 and val <= 360) or val >= 0 and val < 11.25): return "N"
# 1/2 degree increment between the directions
i = 11.25;
for dir in dirs:
if val >= i and val < (i + 22.5):
return dir
i += 22.5
if __name__ == '__main__':
import doctest
doctest.testmod()
from datetime import datetime
from pycdf import CDF
from numpy import diff, array, concatenate
from metobs import mytime
from metobs import buoy
nc_file_pats = dict(buoy='/beach/Instrument_Data/METOBS/MENDOTA/Buoy/%Y/%m/mendota_buoy.%Y-%m-%d.nc',
tower='/beach/Instrument_Data/METOBS/RIG/Tower/%Y/%m/rig_tower.%Y-%m-%d.nc')
def get_fp(inst, name, start, end):
out_file_pat = '%(inst)s_%(var)s_%(start)s_%(end)s.deltas'
fn = out_file_pat % dict(inst=inst,
var=name,
start=start.strftime('%Y-%m-%d'),
end=end.strftime('%Y-%m-%d'))
return open(fn, 'at')
def gen_deltas(inst, start, end, vars):
"""Create a generator for deltas for record variables.
:param inst: Name of the instrument from `nc_file_pats`
:param start: Starting date, inclusive.
:param end: Ending date, exclusive.
:param vars: Variable names from the netcdf files. If None all record
variables are assumed.
"""
dates = mytime.gen_days(start, end)
for dt in dates:
nc = CDF(datetime.strftime(dt, nc_file_pats[inst]))
if not vars:
vars = nc.variables().keys()
my_deltas = {}
for name in nc.variables().keys():
var = nc.var(name)
if var.isrecord() and not name.startswith('qc'):
data = var[:]
my_deltas[name] = diff(data)
yield dt, my_deltas
def to_files(inst, start, end, vars):
for dt, my_deltas in gen_deltas(inst, start, end, vars):
print "writing deltas for", dt
for name in my_deltas:
fp = get_fp(inst, name, start, end)
fp.write('\n'.join([str(v) for v in my_deltas[name]]) + '\n')
fp.close()
if __name__ == '__main__':
from optparse import OptionParser
parser = OptionParser(usage="""%prog [options] (buoy|tower) YYYY-MM-DD YYYY-MM-DD [<var> <var> ...])
If no variables are provided ALL variables will be assumed.
""")
options, args = parser.parse_args()
fmt = '%Y-%m-%d'
try:
inst = args[0]
assert inst in nc_file_pats, "Invalid instrument: %s" % inst
start = mytime.set_tz(datetime.strptime(args[1], fmt))
end = mytime.set_tz(datetime.strptime(args[2], fmt))
if len(args) == 3:
vars = []
else:
vars = args[3:]
except Exception, e:
print e
parser.print_usage()
import sys
sys.exit(1)
to_files(inst, start, end, vars)
"""Solar position calculator and routines.
Ported from Java code which was based on NOAA code at:
http://www.srrb.noaa.gov/highlights/sunrise/azel.html
This was specifically ported to retrieve values for the maximum solar flux::
(cosZen * solarConstant)
Simple usage is as follows::
>>> from metobs.data import solar_position as sp
>>> from datetime import datetime
>>> date = datetime.now()
# or date = datetime(2009,4,1,12,0,0)
>>> s = sp.solar_position(date, 45, -90)
>>> print s
<SolarPosition
date: 2009-04-01 14:11:28.931523-05:00
latitude: 45
longitude: -90
earthSunDistance: 0.999456181732
azimuth: 205.27957233
solarZen: 42.8125188068
cosZen: 0.733581392872
elevations: 47.1874811932
solarDec: 4.81517855635
eqTime: -3.71773858624
solarConstant: 1369.49188965>
The 'SolarPosition' object is a list like structure that allows you to
access the data in a list like manner::
>>> for val in s:
print val
...:
45
-90
0.999456181732
205.27957233
42.8125188068
0.733581392872
47.1874811932
4.81517855635
-3.71773858624
1369.49188965
>>> s[1]
45
>>> s[sp.SP_COS_ZEN]
0.73358139287238222
Additionally, there are property attributes for each value::
>>> s.solarDec
4.8151785563547014
>>> s.earthSunDistance
0.99945618173240092
Doctests are included to verify all output values are within _ACCEPTABLE_DELTA of
the values generated by the Java version.
https://cvs.ssec.wisc.edu/cgi-bin/cvsweb.cgi/java/tower/SolarPosition.java
"""
__docformat__ = 'restructuredtext en'
__version__ = '$Revision: 1.6 $'
# $Source: /cvsroot/TOOLS/dev/metobs/python/data/metobs/data/solar_position.py,v $
from metobs import mytime
from metobs.data import NaN, approx_eq
import math
_days = ["Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday"]
_months = ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"]
_ACCEPTABLE_DELTA = .0000001
def calcDayOfYear(m, d, leap):
"""Find numerical day-of-year from month, day and leap year.
"""
k = leap and 1 or 2
doy = math.floor((275 * m) / 9) - k * math.floor((m + 9) / 12) + d - 30
return doy
def calcDayOfWeek(juld):
"""
>>> jd = 2440587.5
>>> calcDayOfWeek(jd)
'Thursday'
"""
A = (juld + 1.5) % 7
return _days[int(A)]
def calcJD(year, month, day):
"""Get the interval of time in days and fractions of a day, since January 1,
4713 BC Greenwich noon, Julian proleptic calendar.
>>> calcJD(1970,1,1)
2440587.5
"""
if month <= 2:
year -= 1
month += 12
A = math.floor(year/100)
B = 2 - A + math.floor(A/4)
return (math.floor(365.25 * (year + 4716))
+ math.floor(30.6001 * (month + 1))
+ day + B - 1524.5)
def calcDateFromJD(jd):
"""
>>> jd = 2440587.5
>>> calcDateFromJD(jd)
'1-Jan-1970'
"""
z = math.floor(jd + 0.5)
f = (jd + 0.5) - z
if z < 2299161:
A = z
else:
alpha = math.floor((z - 1867216.25)/36524.25)
A = z + 1 + alpha - math.floor(alpha / 4)
B = A + 1524
C = math.floor((B - 122.1) / 365.25)
D = math.floor(365.25 * C)
E = math.floor((B - D) / 30.6001)
day = int(B - D - math.floor(30.6001 * E) + f)
month = int((E < 14) and E - 1 or E - 13)
year = int((month > 2) and C - 4716 or C - 4715)
return "%d-%s-%d" % (day, _months[month-1], year)
def calcDayFromJD(jd):
"""
>>> jd = 2440587.5
>>> calcDayFromJD(jd)
'1-Jan'
"""
s = calcDateFromJD(jd)
return s[:s.rfind('-')]
def calcTimeJulianCent(jd):
"""
>>> jd = 2440587.5
>>> assert approx_eq(calcTimeJulianCent(jd), -.3, _ACCEPTABLE_DELTA)
"""
return (jd - 2451545.0) / 36525.0
def calcJDFromJulianCent(t):
"""
>>> calcJDFromJulianCent(.3)
2462502.5
"""
return t * 36525.0 + 2451545.0
def calcGeomMeanLongSun(t):
"""
>>> assert approx_eq(calcGeomMeanLongSun(.3), 280.69743628799733, _ACCEPTABLE_DELTA)
"""
L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t)
while L0 > 360.0:
L0 -= 360.0
while L0 < 0.0:
L0 += 360.0
return L0
def calcGeomMeanAnomalySun(t):
"""
>>> assert approx_eq(calcGeomMeanAnomalySun(.3), 11157.244183166998, _ACCEPTABLE_DELTA)
"""
return 357.52911 + t * (35999.05029 - 0.0001537 * t)
def calcEccentricityEarthOrbit(t):
"""
>>> assert approx_eq(calcEccentricityEarthOrbit(.3), 0.016696011497, _ACCEPTABLE_DELTA)
"""
return 0.016708634 - t * (0.000042037 + 0.0000001267 * t)
def calcSunEqOfCenter(t):
"""
>>> assert approx_eq(calcSunEqOfCenter(.3), -0.09394256323343497, _ACCEPTABLE_DELTA)
"""
m = calcGeomMeanAnomalySun(t)
mrad = math.radians(m)
sinm = math.sin(mrad)
sin2m = math.sin(mrad * 2)
sin3m = math.sin(mrad * 3)
return (sinm * (1.914602 - t * (0.004817 + 0.000014 * t))
+ sin2m * (0.019993 - 0.000101 * t)
+ sin3m * 0.000289)
def calcSunTrueLong(t):
"""
>>> assert approx_eq(calcSunTrueLong(.3), 280.6034937247639, _ACCEPTABLE_DELTA)
"""
l0 = calcGeomMeanLongSun(t)
c = calcSunEqOfCenter(t)
return l0 + c
def calcSunTrueAnomaly(t):
"""
>>> assert approx_eq(calcSunTrueAnomaly(.3), 11157.150240603765, _ACCEPTABLE_DELTA)
"""
m = calcGeomMeanAnomalySun(t)
c = calcSunEqOfCenter(t)
return m + c
def calcSunRadVector(t):
"""
>>> assert approx_eq(calcSunRadVector(.3), .9833249591481789, _ACCEPTABLE_DELTA)
"""
v = calcSunTrueAnomaly(t)
e = calcEccentricityEarthOrbit(t)
return (1.000001018 * (1 - e * e)) / (1 + e * math.cos(math.radians(v)))
def calcSunApparentLong(t):
"""
>>> assert approx_eq(calcSunApparentLong(.3), 280.60256404613983, _ACCEPTABLE_DELTA)
"""
o = calcSunTrueLong(t)
omega = 125.04 - 1934.136 * t
return o - 0.00569 - 0.00478 * math.sin(math.radians(omega))
def calcMeanObliquityOfEcliptic(t):
"""
>>> assert approx_eq(calcMeanObliquityOfEcliptic(.3), 23.43538985995861, _ACCEPTABLE_DELTA)
"""
seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)))
return 23.0 + (26.0 + (seconds / 60.0)) / 60.0
def calcObliquityCorrection(t):
"""
>>> assert approx_eq(calcObliquityCorrection(.3), 23.435157804956095, _ACCEPTABLE_DELTA)
"""
e0 = calcMeanObliquityOfEcliptic(t)
omega = 125.04 - 1934.136 * t
return e0 + 0.00256 * math.cos(math.radians(omega))
def calcSunRtAscension(t):
"""
>>> assert approx_eq(calcSunRtAscension(.3), -78.46872825491317, _ACCEPTABLE_DELTA)
"""
e = calcObliquityCorrection(t)
lm = calcSunApparentLong(t)
tananum = (math.cos(math.radians(e)) * math.sin(math.radians(lm)))
tanadenom = (math.cos(math.radians(lm)))
return math.degrees(math.atan2(tananum, tanadenom))
def calcSunDeclination(t):
"""
>>> assert approx_eq(calcSunDeclination(.3), -23.011812303124596, _ACCEPTABLE_DELTA)
"""
e = calcObliquityCorrection(t)
lm = calcSunApparentLong(t)
sint = math.sin(math.radians(e)) * math.sin(math.radians(lm))
return math.degrees(math.asin(sint))
def calcEquationOfTime(t):
"""
>>> assert approx_eq(calcEquationOfTime(.3), -3.3355410113272552, _ACCEPTABLE_DELTA)
"""
epsilon = calcObliquityCorrection(t)
l0 = calcGeomMeanLongSun(t)
e = calcEccentricityEarthOrbit(t)
m = calcGeomMeanAnomalySun(t)
y = math.tan(math.radians(epsilon)/2.0)
y *= y
sin2l0 = math.sin(2.0 * math.radians(l0))
sinm = math.sin(math.radians(m))
cos2l0 = math.cos(2.0 * math.radians(l0))
sin4l0 = math.sin(4.0 * math.radians(l0))
sin2m = math.sin(2.0 * math.radians(m))
Etime = (y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0
- 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m)
return math.degrees(Etime) * 4.0
SP_DATE = 0
SP_LAT = 1
SP_LON = 2
SP_EARTH_SUN_DIST = 3
SP_AZMTH = 4
SP_SLR_ZEN = 5
SP_COS_ZEN = 6
SP_ELEV = 7
SP_SLR_DEC = 8
SP_EQ_TIME = 9
SP_SLR_CONST = 10
SP_MAX_FLUX = 11
class SolarPosition(object):
VAL_CNT = 12
def __init__(self, lat, lon, earthRadVec, azimuth, solarZen, coszen,
elevation, solarDec, eqtime, date):
self._lat = lat
self._lon = lon
self._earthRadVec = earthRadVec
self._azimuth = azimuth
self._solarZen = solarZen
self._coszen = coszen
self._elevation = elevation
self._solarDec = solarDec
self._eqTime = eqtime
self._date = date
self._cur_idx = 0
@property
def latitude(self):
return self._lat
@property
def longitude(self):
return self._lon
@property
def earthSunDistance(self):
return self._earthRadVec
@property
def azimuth(self):
return self._azimuth
@property
def solarZen(self):
return self._solarZen
@property
def cosZen(self):
return self._coszen
@property
def elevation(self):
return self._elevation
@property
def solarDec(self):
return self._solarDec
@property
def eqTime(self):
return self._eqTime
@property
def timenow(self):
return mytime.datetime_to_epoch(self._date) / 3600
@property
def date(self):
return self._date
@property
def solarConstant(self):
return 1368.0 * math.pow(1.000001018/self.earthSunDistance, 2)
@property
def maxFlux(self):
return self.cosZen() * self.solarConstant()
#
# list methods
#
def __len__(self): return self.VAL_CNT
def __getitem__(self, idx):
if idx == 0: return self.date
if idx == 1: return self.latitude
if idx == 2: return self.longitude
if idx == 3: return self.earthSunDistance
if idx == 4: return self.azimuth
if idx == 5: return self.solarZen
if idx == 6: return self.cosZen
if idx == 7: return self.elevation
if idx == 8: return self.solarDec
if idx == 9: return self.eqTime
if idx == 10: return self.solarConstant
if idx == 11: return self.maxFlux
raise IndexError
#
# Iterator implementation
#
def __iter__(self): return self
def next(self):
if self._cur_idx == self.VAL_CNT:
self._cur_idx = 0
raise StopIteration
val = self[self._cur_idx]
self._cur_idx += 1
return val
def __str__(self):
return """<SolarPosition
date: %s
latitude: %s
longitude: %s
earthSunDistance: %s
azimuth: %s
solarZen: %s
cosZen: %s
elevations: %s
solarDec: %s
eqTime: %s
solarConstant: %s
maxFlux: %s>""" % tuple([x for x in self])
def solar_position(date, latitude, longitude):
"""Get solar position parameters.
- `date`: date for which to generate parameters. If 'date' has not tzinfo
it assumed to be in localtime.
- `latitude`: Station latitude
- `longitude`: Station longitude
Test values were results of running the SolarPosition Java code in cvs at
https://cvs.ssec.wisc.edu/cgi-bin/cvsweb.cgi/java/tower/SolarPosition.java
>>> from datetime import datetime
>>> date = datetime(1970,1,1,12,0,0)
>>> s = solar_position(date, 45, -90)
>>> assert approx_eq(s.eqTime, -3.62297385346387, _ACCEPTABLE_DELTA)
>>> assert approx_eq(s.solarDec, -22.99641598015454, _ACCEPTABLE_DELTA)
>>> assert approx_eq(s.azimuth, 179.10076689553773, _ACCEPTABLE_DELTA)
>>> assert approx_eq(s.elevation, 22.038213553352705, _ACCEPTABLE_DELTA)
>>> assert approx_eq(s.cosZen, 0.3752248974802099, _ACCEPTABLE_DELTA)
"""
if not date.tzinfo:
date = mytime.set_tz(date, tz=mytime.LocalTimezone)
year, month, day, hh, mm, ss = date.timetuple()[:6]
tzoffset = date.tzinfo.utcoffset(date)
if tzoffset.days < 0:
GMToffset = -(24 - tzoffset.seconds/3600)
else:
GMToffset = tzoffset.seconds/3600
timenow = hh + mm/60 + ss/3600 - GMToffset
JD = calcJD(year, month, day)
T = calcTimeJulianCent(JD + timenow/24.0)
solarDec = calcSunDeclination(T)
eqTime = calcEquationOfTime(T)
earthRadVec = calcSunRadVector(T)
solarTimeFix = eqTime + 4.0 * longitude - 60.0 * GMToffset
trueSolarTime = hh * 60.0 + mm + ss/60.0 + solarTimeFix
while trueSolarTime > 1440:
trueSolarTime -= 1440
hourAngle = trueSolarTime / 4.0 - 180.0
if hourAngle < -180:
hourAngle += 360.0
haRad = math.radians(hourAngle)
csz = (math.sin(math.radians(latitude)) *
math.sin(math.radians(solarDec)) +
math.cos(math.radians(latitude)) *
math.cos(math.radians(solarDec)) * math.cos(haRad))
if csz > 1.0:
csz = 1.0
elif csz < -1.0:
csz = -1.0
zenith = math.degrees(math.acos(csz))
azDenom = ( math.cos(math.radians(latitude)) *
math.sin(math.radians(zenith)) )
if abs(azDenom) > 0.001:
azRad = ((( math.sin(math.radians(latitude)) *
math.cos(math.radians(zenith)) ) -
math.sin(math.radians(solarDec))) / azDenom)
if abs(azRad) > 1.0:
if azRad < 0:
azRad = -1.0
else:
azRad = 1.0
azimuth = 180.0 - math.degrees(math.acos(azRad))
if hourAngle > 0.0:
azimuth = -azimuth
else:
if latitude > 0.0:
azimuth = 180.0
else:
azimuth = 0.0
if azimuth < 0.0:
azimuth += 360.0
exoatmElevation = 90.0 - zenith
refractionCorrection = 0.0
if (exoatmElevation > 85.0):
refractionCorrection = 0.0
else:
te = math.tan (math.radians(exoatmElevation))
if (exoatmElevation > 5.0):
refractionCorrection = (58.1 / te - 0.07 / (te*te*te)
+ 0.000086 / (te*te*te*te*te))
elif exoatmElevation > -0.575:
refractionCorrection = (1735.0 + exoatmElevation *
(-518.2 + exoatmElevation * (103.4 +
exoatmElevation * (-12.79 +
exoatmElevation * 0.711) ) ))
else:
refractionCorrection = -20.774 / te
refractionCorrection = refractionCorrection / 3600.0
solarZen = zenith - refractionCorrection
elevation = 0.
coszen = 0.
if solarZen < 108.0:
elevation = 90.0 - solarZen
if solarZen < 90.0:
coszen = math.cos(math.radians(solarZen))
else:
azumuth = NaN
elevation = NaN
coszen = NaN
return SolarPosition(latitude, longitude, earthRadVec, azimuth, solarZen, coszen,
elevation, solarDec, eqTime, date)
if __name__ == '__main__':
import doctest
doctest.testmod()
\ No newline at end of file
[egg_info]
tag_build = b
tag_svn_revision = true
\ No newline at end of file
setup.py 0 → 100644
# If true, then the svn revision won't be used to calculate the
# revision (set to True for real releases)
RELEASE = False
from setuptools import setup, find_packages
import sys, os
from datetime import datetime
from glob import glob
from metobs import mytime
classifiers = ""
version = '0.1'
filelist = [x for x in glob('metobs/**/*') if 'CVS' not in x]
def most_recent_modification(filelist=[]):
return max( [datetime.fromtimestamp(os.stat(x).st_mtime, mytime.UTCTimezone()) for x in filelist] ).strftime('%Y%m%d.%H%M')
_requires = ['metobs.mytime', 'fpconst>=0.7.2', 'pydap']
setup(name='metobs.data',
version=version + '.' + most_recent_modification(filelist),
description="",
long_description="""MetObs data processing libraries.""",
classifiers=filter(None, classifiers.split("\n")),
keywords='ssec metobs',
author='Bruce Flynn, SSEC',
author_email='brucef@ssec.wisc.edu',
url='',
packages=find_packages(exclude=['ez_setup', 'examples', 'tests']),
namespace_packages=['metobs'],
include_package_data=True,
zip_safe=False,
install_requires=_requires,
dependency_links = ['http://larch.ssec.wisc.edu/cgi-bin/repos.cgi']
)
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment