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

Last change on this file since 938 was 938, checked in by tnavarro, 12 years ago

added changetime option in parser + ensure new time axis continuity (source of graphic bug for year-long or multiple years files)
File size: 8.4 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
103## Author: AS. uses gauss_kern and blur_image.
104def smooth2d(field, window=10):
105        ## actually blur_image could work with different coeff on x and y
106        if True in np.isnan(field):
107            print "!! ERROR !! Smooth is a disaster with missing values. This will be fixed."
108            exit()
109        if window > 1:  result = blur_image(field,int(window))
110        else:           result = field
111        return result
112
113## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
114def gauss_kern(size, sizey=None):
115        # Returns a normalized 2D gauss kernel array for convolutions
116        size = int(size)
117        if not sizey:
118                sizey = size
119        else:
120                sizey = int(sizey)
121        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
122        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
123        return g / g.sum()
124
125## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
126def blur_image(im, n, ny=None) :
127        # blurs the image by convolving with a gaussian kernel of typical size n.
128        # The optional keyword argument ny allows for a different size in the y direction.
129        g = gauss_kern(n, sizey=ny)
130        improc = sp_signal.convolve(im, g, mode='same')
131        return improc
132
133## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
134def smooth1d(x,window=11,window_type='hanning'):
135    """smooth the data using a window with requested size.
136    This method is based on the convolution of a scaled window with the signal.
137    The signal is prepared by introducing reflected copies of the signal
138    (with the window size) in both ends so that transient parts are minimized
139    in the begining and end part of the output signal.
140    input:
141        x: the input signal
142        window: the dimension of the smoothing window; should be an odd integer
143        window_type: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
144            flat window will produce a moving average smoothing.
145    output:
146        the smoothed signal
147    example:
148    t=linspace(-2,2,0.1)
149    x=sin(t)+randn(len(t))*0.1
150    y=smooth(x)
151    see also:
152    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
153    scipy.signal.lfilter
154    TODO: the window parameter could be the window itself if an array instead of a string   
155    """
156    if True in np.isnan(field):
157        print "!! ERROR !! Smooth is a disaster with missing values. This will be fixed."
158        exit()   
159    x = np.array(x)
160    if x.ndim != 1:
161        raise ValueError, "smooth only accepts 1 dimension arrays."
162    if x.size < window:
163        raise ValueError, "Input vector needs to be bigger than window size."
164    if window<3:
165        return x
166    if not window_type in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
167        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
168    s=np.r_[x[window-1:0:-1],x,x[-1:-window:-1]]
169    #print(len(s))
170    if window == 'flat': #moving average
171        w=np.ones(window,'d')
172    else:
173        w=eval('np.'+window_type+'(window)')
174    y=np.convolve(w/w.sum(),s,mode='valid')
175    return y
176
177########################
178#### TIME CONVERTER ####
179########################
180
181#  mars_sol2ls
182#  author T. Navarro
183#  -----------------
184#  convert a given martian day number (sol)
185#  into corresponding solar longitude, Ls (in degr.),
186#  where sol=0=Ls=0 is the
187#  northern hemisphere spring equinox.
188#  -----------------
189def mars_sol2ls(soltabin,forcecontinuity=True):
190      year_day  = 668.6
191      peri_day  = 485.35
192      e_elips   = 0.09340
193      radtodeg  = 57.2957795130823
194      timeperi  = 1.90258341759902
195      if type(soltabin).__name__ in ['int','float','float32','float64']:
196        soltab=[soltabin]
197        solout=np.zeros([1])
198      else:
199        soltab=soltabin
200        solout=np.zeros([len(soltab)])
201      i=0
202      for sol in soltab:
203         zz=(sol-peri_day)/year_day
204         zanom=2.*np.pi*(zz-np.floor(zz))
205         xref=np.abs(zanom)
206#  The equation zx0 - e * sin (zx0) = xref, solved by Newton
207         zx0=xref+e_elips*m.sin(xref)
208         iter=0
209         while iter <= 10:
210            iter=iter+1
211            zdx=-(zx0-e_elips*m.sin(zx0)-xref)/(1.-e_elips*m.cos(zx0))
212            if(np.abs(zdx) <= (1.e-7)):
213              continue
214            zx0=zx0+zdx
215         zx0=zx0+zdx
216         if(zanom < 0.): zx0=-zx0
217# compute true anomaly zteta, now that eccentric anomaly zx0 is known
218         zteta=2.*m.atan(m.sqrt((1.+e_elips)/(1.-e_elips))*m.tan(zx0/2.))
219# compute Ls
220         ls=zteta-timeperi
221         if(ls < 0.): ls=ls+2.*np.pi
222         if(ls > 2.*np.pi): ls=ls-2.*np.pi
223# convert Ls in deg.
224         ls=radtodeg*ls
225         solout[i]=ls
226         i=i+1
227      if forcecontinuity:
228         for iii in range(len(soltab)-1):
229           while solout[iii+1] - solout[iii] > 180. :  solout[iii+1] = solout[iii+1] - 360.
230           while solout[iii] - solout[iii+1] > 180. :  solout[iii+1] = solout[iii+1] + 360.
231      return solout
232
Note: See TracBrowser for help on using the repository browser.