Commit da2bad53 authored by Miguel Tremblay's avatar Miguel Tremblay
Browse files

Move functions from metro_preprocess_fsint2.py in their relevant places in the metro toolbox.

parent e2b164e4
......@@ -56,6 +56,7 @@ import numpy
import metro_logger
import Sun
from toolbox import metro_physics
from toolbox import metro_constant
from toolbox import metro_util
from toolbox import metro_date
......@@ -109,26 +110,15 @@ class Metro_preprocess_fsint2(Metro_preprocess):
wf_interpolated_data,
station_data):
"""
Name: __set_attribute
Set the attributes needed by this clas
Parameters:[I] metro_data wf_controlled_data : controlled data. Read-only
[I] metro_data wf_interpolated_data : container of the
interpolated data.
Returns: None
Functions Called: self.forecast_data.get_matrix_col
station_data.get_latitude
station_data.get_longitude
__get_eot
Description: Set the attributes needed by this class
Notes: <other information, if any>
Parameters:
wf_controlled_data (metro_data) : controlled data. Read-only
wf_interpolated_data (metro_data) : container of the
interpolated data.
Revision History:
Author Date Reason
Miguel Tremblay July 2nd 2004
"""
Returns: None
"""
# Get the cTime of the beginning of data
nTime = wf_controlled_data.get_matrix_col('FORECAST_TIME')[0]
......@@ -136,7 +126,7 @@ class Metro_preprocess_fsint2(Metro_preprocess):
self.fLat = station_data.get_latitude()
self.fLon = station_data.get_longitude()
# Set the constant for the position of the earth around the sun
(self.fEot, self.fR0r, self.tDeclsc) = self.__get_eot(nTime)
(self.fEot, self.fR0r, self.tDeclsc) = metro_date.get_eot(nTime, self.fLat)
self.__set_sunrise_sunset(wf_controlled_data)
def __set_theoretical_flux(self, wf_controlled_data, \
......@@ -208,8 +198,7 @@ class Metro_preprocess_fsint2(Metro_preprocess):
npAT = wf_controlled_data.get_matrix_col('AT')
npIR = npCoeff1*npAT+npCoeff2
wf_controlled_data.set_matrix_col('IR', npIR)
npIR = metro_util.interpolate(npTime, npIR, \
metro_constant.fTimeStep)
npIR = metro_util.interpolate(npTime, npIR)
wf_interpolated_data.append_matrix_col('IR', npIR)
......@@ -220,11 +209,24 @@ class Metro_preprocess_fsint2(Metro_preprocess):
Parameters:
wf_controlled_data (weather forecast)
"""
# Get data
npTime = wf_controlled_data.get_matrix_col('Time')
npSF = self.__get_sf(wf_controlled_data, self.fSunrise, self.fSunset)
npCloudOctal = wf_controlled_data.get_matrix_col('CC')
npTimeHour = wf_controlled_data.get_matrix_col('Hour')
fStartForecastTime = wf_controlled_data.\
get_matrix_col('FORECAST_TIME')[0]
# Get solar fluxes for this cloud cover for this specific day
npSF = metro_physics.get_sf(npCloudOctal, npTimeHour, \
fStartForecastTime,\
self.fSunrise, self.fSunset,\
self.fLat, self.fLon)
# Set value in matrix
wf_controlled_data.set_matrix_col('SF', npSF)
npSF2 = metro_util.interpolate(npTime, npSF, \
metro_constant.fTimeStep)
# Set value in interpolated matrix.
npSF2 = metro_util.interpolate(npTime, npSF)
wf_interpolated_data.append_matrix_col('SF', npSF2)
......@@ -253,61 +255,6 @@ class Metro_preprocess_fsint2(Metro_preprocess):
self.fSunset = fSunsetTimeUTC
def __get_sf(self, wf_controlled_data, fSunriseTimeUTC, fSunsetTimeUTC):
"""
Description: Return an array containing the values of SF.
Parameters:
wf_controlled_data (weather forecast)
fSunriseTimeUTC (float): sunrise time in UTC
fSunsetTimeUT (float): sunset time in UTC
"""
npTimeHour = wf_controlled_data.get_matrix_col('Hour')
npCloudsOctal = wf_controlled_data.get_matrix_col('CC')
nTimeHourLength = len(npTimeHour)
npSft = numpy.zeros(nTimeHourLength, dtype=numpy.float)
npCoeff = numpy.zeros(nTimeHourLength, dtype=numpy.float)
###### In the night, the solar flux is null ###############
for i in range(0, nTimeHourLength):
# Current hour is needed for the computation of
# fDh in the theoritical solar flux.
nCurrentHour = (npTimeHour[i])%24
# atmospheric forecast is before the sunrise
# or after the sunset
if self.__in_the_dark(nCurrentHour, fSunriseTimeUTC, \
fSunsetTimeUTC):
npSft[i] = 0
else:
# Position of the sun around the earth, in radian
fDh = pi*(nCurrentHour/12.0 + self.fLon/180 - 1) + self.fEot
fCosz = self.tDeclsc[0] + \
self.tDeclsc[1]*cos(fDh)
npSft[i] = max(0.0, fCosz)*self.fR0r
npCoeff = -1.56e-12*npSft**4 + 5.972e-9*npSft**3 -\
8.364e-6*npSft**2 + 5.183e-3*npSft - 0.435
npCoeff = numpy.where(npCoeff > 0, npCoeff, 0.0)
# Set npCloudsPercent to be able to reference it in the
# numpy.where method.
npCloudsPercentDay = npCloudsOctal
# Correction based on the clouds
for i in range(0,9):
nPercentDay = metro_constant.lCloudsDay[i]
npCloudsPercentDay = numpy.where(npCloudsOctal == i,\
nPercentDay, npCloudsPercentDay)
npCloudsPercentDay = numpy.where(npCloudsPercentDay == 0, 1.0, \
npCloudsPercentDay)
# Solar flux
npSF3 = npSft * npCoeff * npCloudsPercentDay
return npSF3
def __get_cloud_coefficient(self, wf_controlled_data):
"""
Get the coefficient D1 and D2 as described in the metro article
......@@ -333,119 +280,5 @@ class Metro_preprocess_fsint2(Metro_preprocess):
return (npCoeff1, npCoeff2)
def __get_eot(self, nTime):
"""
Name: __get_eot
Parameters: [I] int nTime : cTime of the beginning of the forecast
Returns: tuple (double fEot, double fR0r, tuple tDeclsc)
dEot: Correction for the equation of time
dR0r: Corrected solar constant for the equation of time
tDeclsc: Declinaison
Functions Called: __Solcons
time.gmtime
calendar.isleap
cos, sin
Description: Subroutine computing the part of the equation of time
needed in the computing of the theoritical solar flux
Correction originating of the CMC GEM model.
Revision History:
Author Date Reason
Miguel Tremblay June 30th 2004
"""
# Convert ctime to python tuple for time.
# see http://www.python.org/doc/current/lib/module-time.html
tDate = time.gmtime(nTime)
# Julian date is the 7th argument
fJulianDate = tDate[7] + tDate[3]/24.0
# Check if it is a leap year
if(calendar.isleap(tDate[0])):
fDivide = 366.0
else:
fDivide = 365.0
# Correction for "equation of time"
fA = fJulianDate/fDivide*2*pi
self.fR0r = self.__Solcons(fA)*metro_constant.fConsol
fRdecl = 0.412*cos((fJulianDate+10.0)*2.0*pi/fDivide-pi)
fDeclsc1 = sin(self.fLat*pi/180.0)*sin(fRdecl)
fDeclsc2 = cos(self.fLat*pi/180.0)*cos(fRdecl)
self.tDeclsc = (fDeclsc1, fDeclsc2)
# in minutes
self.fEot = 0.002733 -7.343*sin(fA)+ .5519*cos(fA) -9.47*sin(2.0*fA) \
-3.02*cos(2.0*fA) -0.3289*sin(3.*fA) -0.07581*cos(3.0*fA) \
-0.1935*sin(4.0*fA) -0.1245*cos(4.0*fA)
# Express in fraction of hour
self.fEot = self.fEot/60.0
# Express in radians
self.fEot = self.fEot*15*pi/180.0
return (self.fEot, self.fR0r, self.tDeclsc)
def __Solcons(self, dAlf):
"""
Name: __Solcons
Parameters: [I] double dAlf : Solar constant to correct the excentricity
Returns: double dVar : Variation of the solar constant
Functions Called: cos, sin
Description: Statement function that calculates the variation of the
solar constant as a function of the julian day. (dAlf, in radians)
Notes: <other information, if any>
Revision History:
Author Date Reason
Miguel Tremblay June 30th 2004
"""
dVar = 1.0/(1.0-9.464e-4*sin(dAlf)-0.01671*cos(dAlf)- \
+ 1.489e-4*cos(2.0*dAlf)-2.917e-5*sin(3.0*dAlf)- \
+ 3.438e-4*cos(4.0*dAlf))**2
return dVar
def __in_the_dark(self, nCurrentTime, fSunrise, fSunset):
"""
Name: __in_the_dark
Parameters: [I] int nCurrentTime. Current time. In [0,24]
[I] float fSunrise. Sunrise time as returned by Sun.py
Value could be > 24.
[I] float fSunset. Sunset time as returned by Sun.py
Value could be > 24.
Returns: boul bDark
Functions Called: cos, sin
Description: Sometimes, value returned by Sun.py are over 24. Since
the time of day is needed with modulo 24, a special check must be
performed. See https://gna.org/bugs/?8277 for more details.
Notes: <other information, if any>
Revision History:
Author Date Reason
Miguel Tremblay January 15th 2006
"""
bDark=False
if ((fSunset%24 > fSunrise%24) and \
(nCurrentTime < fSunrise or \
nCurrentTime > fSunset)) or \
(not(fSunset%24 > fSunrise%24) and \
(nCurrentTime > fSunset%24 and \
nCurrentTime < fSunrise%24)):
bDark=True
return bDark
......@@ -42,10 +42,18 @@ Notes: Time zone string used through this file correspond to the code
on your system.
"""
import xml.utils.iso8601
import math
from math import pi
from math import sin
from math import cos
# Date import
import time
import datetime
import calendar
import xml.utils.iso8601
import os
import string
......@@ -286,3 +294,94 @@ def str_to_at_least_two_digits(nNumber):
sRes = '0'+sRes
return sRes
def in_the_dark(nCurrentTime, fSunrise, fSunset):
"""
Sometimes, value returned by Sun.py are over 24. Since
the time of day is needed with modulo 24, a special check must be
performed. See https://gna.org/bugs/?8277 for more details.
Parameters:
nCurrentTime (int): Current time. In [0,24]
fSunrise (float): Sunrise time as returned by Sun.py. Value could be > 24.
fSunset (float): Sunset time as returned by Sun.py. Value could be > 24.
Returns bDark (bool): True if between fSunset and fSunrise.
False if between fSunrise and fSunset.
"""
bDark=False
if ((fSunset%24 > fSunrise%24) and \
(nCurrentTime < fSunrise or \
nCurrentTime > fSunset)) or \
(not(fSunset%24 > fSunrise%24) and \
(nCurrentTime > fSunset%24 and \
nCurrentTime < fSunrise%24)):
bDark=True
return bDark
def get_eot(fTime, fLat):
"""
Subroutine computing the part of the equation of time
needed in the computing of the theoritical solar flux
Correction originating of the CMC GEM model.
Parameters:
fTime (float): cTime of the beginning of the forecast
fLat (float): latitude of emplacement
Returns: tuple (double fEot, double fR0r, tuple tDeclsc)
dEot: Correction for the equation of time
dR0r: Corrected solar constant for the equation of time
tDeclsc: Declinaison
"""
# Convert ctime to python tuple for time.
# see http://www.python.org/doc/current/lib/module-time.html
tDate = time.gmtime(fTime)
# Julian date is the 7th argument
fJulianDate = tDate[7] + tDate[3]/24.0
# Check if it is a leap year
if(calendar.isleap(tDate[0])):
fDivide = 366.0
else:
fDivide = 365.0
# Correction for "equation of time"
fA = fJulianDate/fDivide*2*pi
fR0r = solcons(fA)*metro_constant.fConsol
fRdecl = 0.412*cos((fJulianDate+10.0)*2.0*pi/fDivide-pi)
fDeclsc1 = sin(fLat*pi/180.0)*sin(fRdecl)
fDeclsc2 = cos(fLat*pi/180.0)*cos(fRdecl)
tDeclsc = (fDeclsc1, fDeclsc2)
# in minutes
fEot = 0.002733 -7.343*sin(fA)+ .5519*cos(fA) -9.47*sin(2.0*fA) \
-3.02*cos(2.0*fA) -0.3289*sin(3.*fA) -0.07581*cos(3.0*fA) \
-0.1935*sin(4.0*fA) -0.1245*cos(4.0*fA)
# Express in fraction of hour
fEot = fEot/60.0
# Express in radians
fEot = fEot*15*pi/180.0
return (fEot, fR0r, tDeclsc)
def solcons(dAlf):
"""
Statement function that calculates the variation of the
solar constant as a function of the julian day. (dAlf, in radians)
Parameters:
dAlf (double) : Solar constant to correct the excentricity
Returns dVar (double : Variation of the solar constant
"""
dVar = 1.0/(1.0-9.464e-4*sin(dAlf)-0.01671*cos(dAlf)- \
+ 1.489e-4*cos(2.0*dAlf)-2.917e-5*sin(3.0*dAlf)- \
+ 3.438e-4*cos(4.0*dAlf))**2
return dVar
......@@ -42,10 +42,15 @@ Date: August 24th 2004
"""
import math
import numarray
from math import pi
from math import sin
from math import cos
from toolbox import metro_constant
from toolbox import metro_util
import numpy
import metro_constant
import metro_util
import metro_date
def foqst(dTD, dPO):
......@@ -94,3 +99,64 @@ def foew(dPO):
(dPO-35.86+max(0.0,metro_util.sign(28.2,-dDiff))))
return fResult
def get_sf(npCloudsOctal, npTimeHour, fStartForecastTime, \
fSunriseTimeUTC, fSunsetTimeUTC, fLat, fLon):
"""
Description: Return an array containing the values of SF.
Parameters:
npCloudOctal (numpy array): array of octal value representing the
cloud coverage.
fSunriseTimeUTC (float): sunrise time in UTC
fSunsetTimeUT (float): sunset time in UTC
nStartForecastTime (??):
fLat (float): latitude of emplacement
fLon (float): longitude of emplacement
Return npSF (numpy array): Array containing the solar flux.
"""
nTimeHourLength = len(npCloudsOctal)
(fEot, fR0r, tDeclsc) = metro_date.get_eot(fStartForecastTime, fLat)
npSft = numpy.zeros(nTimeHourLength, dtype=numpy.float)
npCoeff = numpy.zeros(nTimeHourLength, dtype=numpy.float)
###### In the night, the solar flux is null ###############
for i in range(0, nTimeHourLength):
# Current hour is needed for the computation of
# fDh in the theoritical solar flux.
nCurrentHour = (npTimeHour[i])%24
# atmospheric forecast is before the sunrise
# or after the sunset
if metro_date.in_the_dark(nCurrentHour, fSunriseTimeUTC, \
fSunsetTimeUTC):
npSft[i] = 0
else:
# Position of the sun around the earth, in radian
fDh = pi*(nCurrentHour/12.0 + fLon/180 - 1) + fEot
fCosz = tDeclsc[0] + tDeclsc[1]*cos(fDh)
npSft[i] = max(0.0, fCosz)*fR0r
npCoeff = -1.56e-12*npSft**4 + 5.972e-9*npSft**3 -\
8.364e-6*npSft**2 + 5.183e-3*npSft - 0.435
npCoeff = numpy.where(npCoeff > 0, npCoeff, 0.0)
# Set npCloudsPercent to be able to reference it in the
# numpy.where method.
npCloudsPercentDay = npCloudsOctal
# Correction based on the clouds
for i in range(0,9):
nPercentDay = metro_constant.lCloudsDay[i]
npCloudsPercentDay = numpy.where(npCloudsOctal == i,\
nPercentDay, npCloudsPercentDay)
npCloudsPercentDay = numpy.where(npCloudsPercentDay == 0, 1.0, \
npCloudsPercentDay)
# Solar flux
npSF = npSft * npCoeff * npCloudsPercentDay
return npSF
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment