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

Big futz around, version bumped to 0.4a

parent 9654c01c
Branches
Tags
No related merge requests found
"""Functions for calculating the dervied data.
All linear conversion support ``numpy.ndarray`` objects.
"""
__version__ = '$Revision: 1.15 $'
__docformat__ = 'restructuredtext en'
import math
from datetime import timedelta
import fpconst
from numpy import array, average, zeros, radians, degrees, arctan2, sin, cos, ndarray, log, minimum, exp, sqrt
from numpy.ma import masked_array, average as masked_average, MaskedArray, masked_where, column_stack
try:
from pydap.client import open_url as dapopen
except:
dapopen = open
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
"""
if type(arr) == ndarray:
avg_func = average
array_func = array
elif type(arr) == MaskedArray:
avg_func = masked_average
array_func = masked_array
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]))
if type(arr) == ndarray: avg_out.fill(NaN)
elif type(arr) == MaskedArray: masked_where(avg_out == 0, avg_out)
# 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:] = avg_func(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:,:]
avg_arr = avg_func(avg_arr, axis=0)
arr_out.append(avg_arr)
if type(arr) == ndarray: arr_out = array_func(arr_out)
elif type(arr) == MaskedArray:
arr_out = column_stack(arr_out).transpose()
return 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 tempC is None or relhum is None:
return NaN
gasconst = 461.5
latheat = 2500800.0
dp = 1.0 / ( 1.0 / ( 273.15 + tempC ) - gasconst * log( (0.0 + relhum) / 100 ) / \
( latheat - tempC * 2397.5 ))
return minimum(dp - 273.15, tempC)
def relhum(airTempK, dewpointTempK):
"""
Algorithm 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 * (pressureMB.max()/pressureMB)**.286
return pT
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 data * float(factor) + offset
def mm2in(val):
"""Convert millimeters to inches.
"""
return linearConvert(val, 0.0393700787)
def in2mm(val):
"""Convert inches to millimeters.
"""
return linearConvert(val, 1/0.0393700787)
def c2f (val):
"""Degrees celsius to fahrenheit.
>>> c2f(0)
32.0
>>> c2f(100)
212.0
"""
return linearConvert(val, (9/5.), 32)
def f2c(val):
"""Degrees fahrenheit to celsius.
>>> f2c(32)
0.0
>>> f2c(212)
100.0
"""
return linearConvert(val-32, (5/9.), 0)
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 knots2mps(val):
"""Speed in knots to meters per second.
"""
return linearConvert(val, 1/1.9438445)
def altimeter(p, alt):
"""Compute altimeter from pressure and altitude.
Converted from code provided by TomW.
:param p: pressure in hPa.
:param alt: altitude of the measurement in meters.
:returns: altimeter in inHg
"""
n = .190284
c1= .0065 * pow(1013.25, n)/288.
c2 = alt / 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
def array_convert(arr, sunits, cunits):
"""Converts the data in array arr from sunits to cunits.
:param arr: Array to be converted
:type arr: ndarray
:param sunits: Starting Units = 'C', 'F', 'm/s', 'knots', 'hpa', 'deg', 'in', 'mm', 'w/m2'
:type sunits: string
:param cunits: Convert to Units = Same choices as sunits
:type cunits: string
"""
if sunits == cunits: return arr
if not isinstance(arr, ndarray): raise ValueError("Array must be of type numpy.ndarray, not %r" % (type(arr),))
if sunits == 'c':
if cunits == 'f':
a = c2f(arr)
return a
elif sunits == 'f':
if cunits == 'c':
a = f2c(arr)
return a
elif sunits == 'm/s':
if cunits == 'knots':
a = mps2knots(arr)
return a
elif sunits == 'knots':
if cunits == 'm/s':
a = knots2mps(arr)
return a
elif sunits == 'in':
if cunits == 'mm':
a = in2mm(arr)
return a
elif sunits == 'mm':
if cunits == 'in':
a = mm2in(arr)
return a
raise ValueError("sunits or cunits was not an acceptable unit")
if __name__ == '__main__':
import doctest
doctest.testmod()
from .time import to_unix_timestamp
from .units import *
from .calc import *
import math
import numpy as np
NaN = float('nan')
is_nan = lambda a: a != a
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 tempC is None or relhum is None:
return NaN
gasconst = 461.5
latheat = 2500800.0
dp = 1.0 / (1.0 / (273.15 + tempC) - gasconst * np.log((0.0 + relhum) / 100) /
(latheat - tempC * 2397.5))
return min(dp - 273.15, tempC)
def relhum(airTempK, dewpointTempK):
"""
Algorithm 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 * (pressureMB.max() / pressureMB) ** .286
return pT
def altimeter(p, alt):
"""Compute altimeter from pressure and altitude.
Converted from code provided by TomW.
:param p: pressure in hPa.
:param alt: altitude of the measurement in meters.
:returns: altimeter in inHg
"""
n = .190284
c1 = .0065 * pow(1013.25, n) / 288.
c2 = alt / 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
def wind_vector_components(windspd, winddir):
"""Decompose scalar or list/array polar wind direction and speed data
into the horizontal and vertical vector components and speed vector.
Inputs can be scalar or arrays.
"""
dir_rad = np.deg2rad(winddir)
spd_arr = np.array(windspd)
V_e = spd_arr * np.sin(dir_rad)
V_n = spd_arr * np.cos(dir_rad)
U_spd = np.sqrt(pow(V_e, 2) + pow(V_n, 2))
return V_e, V_n, U_spd
def wind_vector_degrees(vector_east, vector_north):
"""Re-compose horizontal (east/west) and vertical (north/south) vector
components into wind direction in degrees.
Inputs can be scalar or arrays.
"""
rads = np.arctan2(vector_east, vector_north)
winddir = np.rad2deg(rads)
if isinstance(winddir, np.ndarray):
winddir[np.less(winddir, 0)] += 360
elif winddir < 0:
winddir += 360
return winddir % 360
def mean_wind_vector(windspd, winddir):
V_e, V_n, V_spd = wind_vector_components(windspd, winddir)
avg_dir = wind_vector_degrees(np.mean(V_e), np.mean(V_n))
return avg_dir, np.mean(V_spd)
from datetime import datetime
def test_to_unix_timestamp():
from metobs.data.time import to_unix_timestamp
assert to_unix_timestamp(datetime(1970, 1, 1)) == 0
def test_hhmm_to_secs():
from metobs.data.time import hhmm_to_offset
assert hhmm_to_offset('2400') == 86400, "Can't handle > 23:59"
assert hhmm_to_offset('2401') == 86460, "Can't handle > 23:59"
assert hhmm_to_offset('0') == 0, "Can't handle short times"
assert hhmm_to_offset('001') == 60, "Can't handle leading 0"
import unittest
class MeanWindVectorTests(unittest.TestCase):
def _fut(self, winddir, windspd=None):
from metobs.data.calc import mean_wind_vector
windspd = windspd or [1]*len(winddir)
return mean_wind_vector(windspd, winddir)[0]
def test_spanning_0_degrees(self):
winddir = self._fut([315, 45])
self.assertAlmostEqual(winddir, 0)
def test_spanning_cardinal_directions(self):
self.assertAlmostEqual(self._fut([45, 135]), 90)
self.assertAlmostEqual(self._fut([135, 225]), 180)
self.assertAlmostEqual(self._fut([225, 315]), 270)
self.assertAlmostEqual(self._fut([315, 45]), 0)
def test_all_zeros(self):
self.assertAlmostEqual(self._fut([0, 0]), 0)
def test_zero_windspd(self):
self.assertAlmostEqual(self._fut([0, 0], windspd=[0, 0]), 0)
def test_45s(self):
self.assertAlmostEqual(self._fut([0, 90]), 45)
self.assertAlmostEqual(self._fut([90, 180]), 135)
self.assertAlmostEqual(self._fut([180, 270]), 225)
self.assertAlmostEqual(self._fut([270, 0]), 315)
from datetime import timedelta
from calendar import timegm
def to_unix_timestamp(dtval):
"""Convert a datetime to a unix timestamp.
"""
return timegm(dtval.utctimetuple())
def hhmm_to_offset(hhmm):
"""Convert a string time, possibly with missing hours and minutes, to an
offset of seconds.
"""
hhmm = '{:04d}'.format(int(hhmm))
return timedelta(hours=int(hhmm[0:2]),
minutes=int(hhmm[2:])).total_seconds()
from .util import linear_convert
def mm2in(val):
"""Convert millimeters to inches.
"""
return linear_convert(val, 0.0393700787)
def in2mm(val):
"""Convert inches to millimeters.
"""
return linear_convert(val, 1 / 0.0393700787)
def c2f(val):
"""Degrees celsius to fahrenheit.
>>> c2f(0)
32.0
>>> c2f(100)
212.0
"""
return linear_convert(val, (9 / 5.), 32)
def f2c(val):
"""Degrees fahrenheit to celsius.
>>> f2c(32)
0.0
>>> f2c(212)
100.0
"""
return linear_convert(val - 32, (5 / 9.), 0)
def mps2mph(val):
"""Speed in meters per second to miles per hour.
"""
return linear_convert(val, 2.23694)
def mps2knots(val):
"""Meters per second to knots.
"""
return linear_convert(val, 1.9438445)
def knots2mps(val):
"""Speed in knots to meters per second.
"""
return linear_convert(val, 1 / 1.9438445)
\ No newline at end of file
from datetime import datetime, timedelta
from .units import *
from .time import to_unix_timestamp
import numpy as np
from numpy.ma import (masked_array, average as masked_average, MaskedArray,
masked_where, column_stack)
NaN = float('nan')
def linear_convert(data, factor=1, offset=0):
"""
>>> linear_convert(0)
0.0
>>> linear_convert(1)
1.0
>>> linear_convert(1, 1, 0)
1.0
>>> linear_convert(1, 90, 10)
100.0
"""
return data * float(factor) + offset
def array_convert(arr, sunits, cunits):
"""Converts the data in array arr from sunits to cunits.
:param arr: Array to be converted
:type arr: ndarray
:param sunits: Starting Units = 'C', 'F', 'm/s', 'knots', 'hpa', 'deg', 'in', 'mm', 'w/m2'
:param cunits: Convert to Units = Same choices as sunits
"""
if sunits == cunits: return arr
if not isinstance(arr, np.ndarray):
raise ValueError("Array must be of type numpy.ndarray, not %r" % (type(arr),))
if sunits == 'c':
if cunits == 'f':
a = c2f(arr)
return a
elif sunits == 'f':
if cunits == 'c':
a = f2c(arr)
return a
elif sunits == 'm/s':
if cunits == 'knots':
a = mps2knots(arr)
return a
elif sunits == 'knots':
if cunits == 'm/s':
a = knots2mps(arr)
return a
elif sunits == 'in':
if cunits == 'mm':
a = in2mm(arr)
return a
elif sunits == 'mm':
if cunits == 'in':
a = mm2in(arr)
return a
raise ValueError("sunits or cunits was not an acceptable unit")
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
"""
if isinstance(arr, np.ndarray):
avg_func = np.average
array_func = np.array
elif isinstance(arr, MaskedArray):
avg_func = masked_average
array_func = masked_array
else:
raise ValueError()
arr_out = []
stop_dt = basetime + timedelta(seconds=interval)
stop_idx = 1
start_idx = 0
# advance through data until there is a timestamp in the averaging interval
while datetime.utcfromtimestamp(arr[0][0]) >= stop_dt:
stop_dt += timedelta(seconds=interval)
for row_idx in range(arr.shape[0]):
row = arr[row_idx]
cur_dt = datetime.utcfromtimestamp(row[0])
# determine stop index
if cur_dt < stop_dt:
stop_idx += 1
# stop and average
else:
# matrix of data to be averaged
avg_input = arr[start_idx:stop_idx, 1:]
avg_out = np.array(np.zeros(arr[0].shape[0]))
if isinstance(arr, np.ndarray):
avg_out.fill(NaN)
elif isinstance(arr, MaskedArray):
masked_where(avg_out == 0, avg_out)
# timestamp for center of averaging interval
int_dt = stop_dt - timedelta(seconds=(interval / 2))
avg_out[0] = to_unix_timestamp(int_dt)
avg_out[1:] = avg_func(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:, :]
avg_arr = avg_func(avg_arr, axis=0)
arr_out.append(avg_arr)
if isinstance(arr, np.ndarray):
arr_out = array_func(arr_out)
elif type(arr) == MaskedArray:
arr_out = column_stack(arr_out).transpose()
return arr_out
[egg_info]
tag_build = b
tag_svn_revision = true
\ No newline at end of file
# 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
classifiers = ""
version = '0.3'
filelist = [x for x in glob('metobs/**/*') if 'CVS' not in x]
_requires = ['metobs.mytime', 'fpconst>=0.7.2']
setup(name='metobs.data',
version=version,
version='0.4a',
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']),
url="",
packages=find_packages(exclude=["*.tests", "*.tests.*", "tests.*", "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']
install_requires=[
'numpy'
],
dependency_links=['http://larch.ssec.wisc.edu/cgi-bin/repos.cgi']
)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment