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

Last change on this file since 1004 was 1002, checked in by aslmd, 12 years ago

UTIL PYTHON planetoplot_v2. possibility to load a field in a restricted interval with e.g. -t 0,1 without performing computations (self.compute=nothing). use for histograms, see histo.py. added smooth2diter. bug correction: ensure numpy array, vector, minmax equal.
File size: 10.0 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################
99#### SMOOTH ####
100################
101### TBD: works with missing values
102
103def smooth2diter(field,n=1):
104        count = 0
105        result = field
106        while count < n:
107            result = smooth2d(result, window=2)
108            count = count + 1
109        return result
110
111## Author: AS. uses gauss_kern and blur_image.
112def smooth2d(field, window=10):
113        ## actually blur_image could work with different coeff on x and y
114        if True in np.isnan(field):
115            print "!! ERROR !! Smooth is a disaster with missing values. This will be fixed."
116            exit()
117        if window > 1:  result = blur_image(field,int(window))
118        else:           result = field
119        return result
120
121## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
122def gauss_kern(size, sizey=None):
123        # Returns a normalized 2D gauss kernel array for convolutions
124        size = int(size)
125        if not sizey:
126                sizey = size
127        else:
128                sizey = int(sizey)
129        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
130        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
131        return g / g.sum()
132
133## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
134def blur_image(im, n, ny=None) :
135        # blurs the image by convolving with a gaussian kernel of typical size n.
136        # The optional keyword argument ny allows for a different size in the y direction.
137        g = gauss_kern(n, sizey=ny)
138        improc = sp_signal.convolve(im, g, mode='same')
139        return improc
140
141## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
142def smooth1d(x,window=11,window_type='hanning'):
143    """smooth the data using a window with requested size.
144    This method is based on the convolution of a scaled window with the signal.
145    The signal is prepared by introducing reflected copies of the signal
146    (with the window size) in both ends so that transient parts are minimized
147    in the begining and end part of the output signal.
148    input:
149        x: the input signal
150        window: the dimension of the smoothing window; should be an odd integer
151        window_type: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
152            flat window will produce a moving average smoothing.
153    output:
154        the smoothed signal
155    example:
156    t=linspace(-2,2,0.1)
157    x=sin(t)+randn(len(t))*0.1
158    y=smooth(x)
159    see also:
160    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
161    scipy.signal.lfilter
162    TODO: the window parameter could be the window itself if an array instead of a string   
163    """
164    if True in np.isnan(field):
165        print "!! ERROR !! Smooth is a disaster with missing values. This will be fixed."
166        exit()   
167    x = np.array(x)
168    if x.ndim != 1:
169        raise ValueError, "smooth only accepts 1 dimension arrays."
170    if x.size < window:
171        raise ValueError, "Input vector needs to be bigger than window size."
172    if window<3:
173        return x
174    if not window_type in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
175        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
176    s=np.r_[x[window-1:0:-1],x,x[-1:-window:-1]]
177    #print(len(s))
178    if window == 'flat': #moving average
179        w=np.ones(window,'d')
180    else:
181        w=eval('np.'+window_type+'(window)')
182    y=np.convolve(w/w.sum(),s,mode='valid')
183    return y
184
185########################
186#### TIME CONVERTER ####
187########################
188
189#  mars_sol2ls
190#  author T. Navarro
191#  -----------------
192#  convert a given martian day number (sol)
193#  into corresponding solar longitude, Ls (in degr.),
194#  where sol=0=Ls=0 is the
195#  northern hemisphere spring equinox.
196#  -----------------
197def mars_sol2ls(soltabin,forcecontinuity=True):
198      year_day  = 668.6
199      peri_day  = 485.35
200      e_elips   = 0.09340
201      radtodeg  = 57.2957795130823
202      timeperi  = 1.90258341759902
203      if type(soltabin).__name__ in ['int','float','float32','float64']:
204        soltab=[soltabin]
205        solout=np.zeros([1])
206      else:
207        soltab=soltabin
208        solout=np.zeros([len(soltab)])
209      i=0
210      for sol in soltab:
211         zz=(sol-peri_day)/year_day
212         zanom=2.*np.pi*(zz-np.floor(zz))
213         xref=np.abs(zanom)
214#  The equation zx0 - e * sin (zx0) = xref, solved by Newton
215         zx0=xref+e_elips*m.sin(xref)
216         iter=0
217         while iter <= 10:
218            iter=iter+1
219            zdx=-(zx0-e_elips*m.sin(zx0)-xref)/(1.-e_elips*m.cos(zx0))
220            if(np.abs(zdx) <= (1.e-7)):
221              continue
222            zx0=zx0+zdx
223         zx0=zx0+zdx
224         if(zanom < 0.): zx0=-zx0
225# compute true anomaly zteta, now that eccentric anomaly zx0 is known
226         zteta=2.*m.atan(m.sqrt((1.+e_elips)/(1.-e_elips))*m.tan(zx0/2.))
227# compute Ls
228         ls=zteta-timeperi
229         if(ls < 0.): ls=ls+2.*np.pi
230         if(ls > 2.*np.pi): ls=ls-2.*np.pi
231# convert Ls in deg.
232         ls=radtodeg*ls
233         solout[i]=ls
234         i=i+1
235      if forcecontinuity:
236         for iii in range(len(soltab)-1):
237           while solout[iii+1] - solout[iii] > 180. :  solout[iii+1] = solout[iii+1] - 360.
238           while solout[iii] - solout[iii+1] > 180. :  solout[iii+1] = solout[iii+1] + 360.
239      return solout
240
241# mars_date
242# author A. Spiga
243# ------------------------------
244# get Ls, sol, utc from a string with format yyyy-mm-dd_hh:00:00
245# -- argument timechar is a vector of such strings indicating dates
246# -- example: timechar = nc.variables['Times'][:] in mesoscale files
247# NB: uses mars_sol2ls function above
248# ------------------------------
249def mars_date(timechar):
250    # some preliminary information
251    days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56]
252    plus_in_month = [ 0, 61,127,193,258,318,372,422,468,515,562,613]
253    # get utc and sol from strings
254    utc = [] ; sol = []
255    for zetime in timechar:
256        dautc = float(zetime[11]+zetime[12]) + float(zetime[14]+zetime[15])/37.
257        dasol = dautc / 24.
258        dasol = dasol + plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9])
259        dasol = dasol - 1 ##les sols GCM commencent a 0
260        utc.append(dautc)
261        sol.append(dasol)
262    sol = np.array(sol)
263    utc = np.array(utc)
264    # get ls from sol
265    ls = mars_sol2ls(sol)
266    return ls, sol, utc
267
268# timecorrect
269# author A. Spiga
270# -----------------------------
271# ensure time axis is monotonic
272# correct negative values
273# -----------------------------
274def timecorrect(time):
275    for ind in range(len(time)-1):
276        if time[ind] < 0.: time[ind] = time[ind] + 24.
277        if time[ind+1] < time[ind]: time[ind+1] = time[ind+1] + 24.
278    return time
279
Note: See TracBrowser for help on using the repository browser.