source: trunk/UTIL/PYTHON/planetoplot_v2/ppcompute.py @ 1176

Last change on this file since 1176 was 1007, checked in by aslmd, 11 years ago

UTIL PYTHON. various updates on example + mcd + plot scripts. nothing major.
File size: 10.7 KB
Line 
1###############################################
2## PLANETOPLOT                               ##
3## --> PPCOMPUTE                             ##
4###############################################
5# python built-in librairies
6import os
7# added librairies
8import numpy as np
9import math  as m
10import scipy.signal as sp_signal
11###############################################
12
13## first a useful function to find settings in a folder in PYTHONPATH
14def findset(whereset,string="planetoplot_v2"):
15    # ... set a default whereset if it was set to None
16    # ... default is in the planetoplot_v2 folder
17    if whereset is None:
18        for path in os.environ['PYTHONPATH'].split(os.pathsep):
19            if string in path: whereset = path
20        if whereset is None: 
21            print "!! ERROR !! "+ string + "not in $PYTHONPATH"
22            print "--> either put it in $PYTHONPATH or change whereset"
23            exit()
24    # ... if the last / is missing put it
25    if whereset[-1] != "/": whereset = whereset + "/"
26    return whereset
27
28##########################
29#### MAX MEAN MIN SUM ####
30#####################################
31#### WITH SUPPORT FOR NaN VALUES ####
32#####################################
33
34## compute min
35## author AS + AC
36def min (field,axis=None): 
37        if field is None: return None
38        if type(field).__name__=='MaskedArray':
39              field.set_fill_value(np.NaN)
40              return np.ma.array(field).min(axis=axis)
41        elif (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
42              return np.ma.masked_invalid(field).min(axis=axis)
43        else: return np.array(field).min(axis=axis)
44
45## compute max
46## author AS + AC
47def max (field,axis=None):
48        if field is None: return None
49        if type(field).__name__=='MaskedArray':
50              field.set_fill_value(np.NaN)
51              return np.ma.array(field).max(axis=axis)
52        elif (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
53              return np.ma.masked_invalid(field).max(axis=axis) 
54        else: return np.array(field).max(axis=axis)
55
56## compute mean
57## author AS + AC
58def mean (field,axis=None):
59        if field is None: return None
60        else: 
61           if type(field).__name__=='MaskedArray':
62              field.set_fill_value(np.NaN)
63              zout=np.ma.array(field).mean(axis=axis)
64              if axis is not None:
65                 zout.set_fill_value(np.NaN)
66                 return zout.filled()
67              else:return zout
68           elif (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
69              zout=np.ma.masked_invalid(field).mean(axis=axis)
70              if axis is not None:
71                 zout.set_fill_value([np.NaN])
72                 return zout.filled()
73              else:return zout
74           else: 
75              return np.array(field).mean(axis=axis)
76
77## compute sum
78## author AS + AC
79def sum (field,axis=None):
80        if field is None: return None
81        else:
82           if type(field).__name__=='MaskedArray':
83              field.set_fill_value(np.NaN)
84              zout=np.ma.array(field).sum(axis=axis)
85              if axis is not None:
86                 zout.set_fill_value(np.NaN)
87                 return zout.filled()
88              else:return zout
89           elif (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
90              zout=np.ma.masked_invalid(field).sum(axis=axis)
91              if axis is not None:
92                 zout.set_fill_value([np.NaN])
93                 return zout.filled()
94              else:return zout
95           else:
96              return np.array(field).sum(axis=axis)
97
98## compute mean over bins
99## author AS
100def meanbin(y,x,bins):
101    bins.sort() # sorting is needed for binning
102    meanvalues = []
103    for iii in range(len(bins)):
104        ## GET VALUES OVER RELEVANT INTERVALS
105        if iii == 0:
106          ind = x<bins[iii+1] ; ys = y[ind]
107        elif iii == len(bins)-1:
108          ind = x>=bins[iii-1] ; ys = y[ind]
109        else:
110          ind = x>=bins[iii-1] ; interm = x[ind] ; intermf = y[ind]
111          ind = interm<bins[iii+1] ; xs = interm[ind] ; ys = intermf[ind]
112        ## COMPUTE MEAN and TREAT NAN CASE
113        meanvalues.append(mean(ys))
114    ## RETURN A NUMPY ARRAY
115    meanvalues = np.array(meanvalues)
116    return meanvalues
117
118################
119#### SMOOTH ####
120################
121### TBD: works with missing values
122
123def smooth2diter(field,n=1):
124        count = 0
125        result = field
126        while count < n:
127            result = smooth2d(result, window=2)
128            count = count + 1
129        return result
130
131## Author: AS. uses gauss_kern and blur_image.
132def smooth2d(field, window=10):
133        ## actually blur_image could work with different coeff on x and y
134        if True in np.isnan(field):
135            print "!! ERROR !! Smooth is a disaster with missing values. This will be fixed."
136            exit()
137        if window > 1:  result = blur_image(field,int(window))
138        else:           result = field
139        return result
140
141## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
142def gauss_kern(size, sizey=None):
143        # Returns a normalized 2D gauss kernel array for convolutions
144        size = int(size)
145        if not sizey:
146                sizey = size
147        else:
148                sizey = int(sizey)
149        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
150        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
151        return g / g.sum()
152
153## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
154def blur_image(im, n, ny=None) :
155        # blurs the image by convolving with a gaussian kernel of typical size n.
156        # The optional keyword argument ny allows for a different size in the y direction.
157        g = gauss_kern(n, sizey=ny)
158        improc = sp_signal.convolve(im, g, mode='same')
159        return improc
160
161## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
162def smooth1d(x,window=11,window_type='hanning'):
163    """smooth the data using a window with requested size.
164    This method is based on the convolution of a scaled window with the signal.
165    The signal is prepared by introducing reflected copies of the signal
166    (with the window size) in both ends so that transient parts are minimized
167    in the begining and end part of the output signal.
168    input:
169        x: the input signal
170        window: the dimension of the smoothing window; should be an odd integer
171        window_type: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
172            flat window will produce a moving average smoothing.
173    output:
174        the smoothed signal
175    example:
176    t=linspace(-2,2,0.1)
177    x=sin(t)+randn(len(t))*0.1
178    y=smooth(x)
179    see also:
180    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
181    scipy.signal.lfilter
182    TODO: the window parameter could be the window itself if an array instead of a string   
183    """
184    if True in np.isnan(field):
185        print "!! ERROR !! Smooth is a disaster with missing values. This will be fixed."
186        exit()   
187    x = np.array(x)
188    if x.ndim != 1:
189        raise ValueError, "smooth only accepts 1 dimension arrays."
190    if x.size < window:
191        raise ValueError, "Input vector needs to be bigger than window size."
192    if window<3:
193        return x
194    if not window_type in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
195        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
196    s=np.r_[x[window-1:0:-1],x,x[-1:-window:-1]]
197    #print(len(s))
198    if window == 'flat': #moving average
199        w=np.ones(window,'d')
200    else:
201        w=eval('np.'+window_type+'(window)')
202    y=np.convolve(w/w.sum(),s,mode='valid')
203    return y
204
205########################
206#### TIME CONVERTER ####
207########################
208
209#  mars_sol2ls
210#  author T. Navarro
211#  -----------------
212#  convert a given martian day number (sol)
213#  into corresponding solar longitude, Ls (in degr.),
214#  where sol=0=Ls=0 is the
215#  northern hemisphere spring equinox.
216#  -----------------
217def mars_sol2ls(soltabin,forcecontinuity=True):
218      year_day  = 668.6
219      peri_day  = 485.35
220      e_elips   = 0.09340
221      radtodeg  = 57.2957795130823
222      timeperi  = 1.90258341759902
223      if type(soltabin).__name__ in ['int','float','float32','float64']:
224        soltab=[soltabin]
225        solout=np.zeros([1])
226      else:
227        soltab=soltabin
228        solout=np.zeros([len(soltab)])
229      i=0
230      for sol in soltab:
231         zz=(sol-peri_day)/year_day
232         zanom=2.*np.pi*(zz-np.floor(zz))
233         xref=np.abs(zanom)
234#  The equation zx0 - e * sin (zx0) = xref, solved by Newton
235         zx0=xref+e_elips*m.sin(xref)
236         iter=0
237         while iter <= 10:
238            iter=iter+1
239            zdx=-(zx0-e_elips*m.sin(zx0)-xref)/(1.-e_elips*m.cos(zx0))
240            if(np.abs(zdx) <= (1.e-7)):
241              continue
242            zx0=zx0+zdx
243         zx0=zx0+zdx
244         if(zanom < 0.): zx0=-zx0
245# compute true anomaly zteta, now that eccentric anomaly zx0 is known
246         zteta=2.*m.atan(m.sqrt((1.+e_elips)/(1.-e_elips))*m.tan(zx0/2.))
247# compute Ls
248         ls=zteta-timeperi
249         if(ls < 0.): ls=ls+2.*np.pi
250         if(ls > 2.*np.pi): ls=ls-2.*np.pi
251# convert Ls in deg.
252         ls=radtodeg*ls
253         solout[i]=ls
254         i=i+1
255      if forcecontinuity:
256         for iii in range(len(soltab)-1):
257           while solout[iii+1] - solout[iii] > 180. :  solout[iii+1] = solout[iii+1] - 360.
258           while solout[iii] - solout[iii+1] > 180. :  solout[iii+1] = solout[iii+1] + 360.
259      return solout
260
261# mars_date
262# author A. Spiga
263# ------------------------------
264# get Ls, sol, utc from a string with format yyyy-mm-dd_hh:00:00
265# -- argument timechar is a vector of such strings indicating dates
266# -- example: timechar = nc.variables['Times'][:] in mesoscale files
267# NB: uses mars_sol2ls function above
268# ------------------------------
269def mars_date(timechar):
270    # some preliminary information
271    days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56]
272    plus_in_month = [ 0, 61,127,193,258,318,372,422,468,515,562,613]
273    # get utc and sol from strings
274    utc = [] ; sol = []
275    for zetime in timechar:
276        dautc = float(zetime[11]+zetime[12]) + float(zetime[14]+zetime[15])/37.
277        dasol = dautc / 24.
278        dasol = dasol + plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9])
279        dasol = dasol - 1 ##les sols GCM commencent a 0
280        utc.append(dautc)
281        sol.append(dasol)
282    sol = np.array(sol)
283    utc = np.array(utc)
284    # get ls from sol
285    ls = mars_sol2ls(sol)
286    return ls, sol, utc
287
288# timecorrect
289# author A. Spiga
290# -----------------------------
291# ensure time axis is monotonic
292# correct negative values
293# -----------------------------
294def timecorrect(time):
295    for ind in range(len(time)-1):
296        if time[ind] < 0.: time[ind] = time[ind] + 24.
297        if time[ind+1] < time[ind]: time[ind+1] = time[ind+1] + 24.
298    return time
299
Note: See TracBrowser for help on using the repository browser.