|
Package spharm ::
Module spharm
|
|
1 """
2 Introduction
3 ============
4
5 This module provides a python interface to the NCAR
6 U{SPHEREPACK<https://www2.cisl.ucar.edu/resources/legacy/spherepack>} library.
7 It is not a one-to-one wrapper for the SPHEREPACK routines, rather
8 it provides a simple interface oriented toward working with
9 atmospheric general circulation model (GCM) data.
10
11 Requirements
12 ============
13 - U{numpy<http://numeric.scipy.org>}, and a fortran compiler
14 supported by numpy.f2py.
15
16 Installation
17 ============
18 - U{Download<http://code.google.com/p/pyspharm/downloads/list>} module source,
19 untar.
20 - run C{python setup.py install} (as root if necessary).
21 The SPHERPACK fortran source files will be downloaded automatically by the
22 setup.py script, since the SPHEREPACK license prohibits redistribution.
23 To specify the fortran compiler to use (e.g. g95) run
24 C{python setup.py config_fc --fcompiler=g95 install}. C{f2py -c --help-fcompiler}
25 will show you what fortran compilers are available.
26
27 Usage
28 =====
29
30 >>> import spharm
31 >>> x=spharm.Spharmt(144,72,rsphere=8e6,gridtype='gaussian',legfunc='computed')
32
33 creates a class instance for spherical harmonic calculations on a 144x72
34 gaussian grid on a sphere with radius 8000 km. The associated legendre
35 functions are recomputed on the fly (instead of pre-computed and stored).
36 Default values of rsphere, gridtype and legfunc are 6.3712e6, 'regular'
37 and 'stored'. Real-world examples are included in the source distribution.
38
39 Class methods
40 =============
41 - grdtospec: grid to spectral transform (spherical harmonic analysis).
42 - spectogrd: spectral to grid transform (spherical harmonic synthesis).
43 - getuv: compute u and v winds from spectral coefficients of vorticity
44 and divergence.
45 - getvrtdivspec: get spectral coefficients of vorticity and divergence
46 from u and v winds.
47 - getgrad: compute the vector gradient given spectral coefficients.
48 - getpsichi: compute streamfunction and velocity potential from winds.
49 - specsmooth: isotropic spectral smoothing.
50
51 Functions
52 =========
53 - regrid: spectral re-gridding, with optional spectral smoothing and/or
54 truncation.
55 - gaussian_lats_wts: compute gaussian latitudes and weights.
56 - getspecindx: compute indices of zonal wavenumber and degree
57 for complex spherical harmonic coefficients.
58 - legendre: compute associated legendre functions.
59 - getgeodesicpts: computes the points on the surface of the sphere
60 corresponding to a twenty-sided (icosahedral) geodesic.
61 - specintrp: spectral interpolation to an arbitrary point on the sphere.
62
63 Conventions
64 ===========
65
66 The gridded data is assumed to be oriented such that i=1 is the
67 Greenwich meridian and j=1 is the northernmost point. Grid indices
68 increase eastward and southward. If nlat is odd the equator is included.
69 If nlat is even the equator will lie half way between points nlat/2
70 and (nlat/2)+1. nlat must be at least 3. For regular grids
71 (gridtype='regular') the poles will be included when nlat is odd.
72 The grid increment in longitude is 2*pi/nlon radians. For example,
73 nlon = 72 for a five degree grid. nlon must be greater than or
74 equal to 4. The efficiency of the computation is improved when nlon
75 is a product of small prime numbers.
76
77 The spectral data is assumed to be in a complex array of dimension
78 (ntrunc+1)*(ntrunc+2)/2. ntrunc is the triangular truncation limit
79 (ntrunc = 42 for T42). ntrunc must be <= nlat-1. Coefficients are
80 ordered so that first (nm=0) is m=0,n=0, second is m=0,n=1,
81 nm=ntrunc is m=0,n=ntrunc, nm=ntrunc+1 is m=1,n=1, etc.
82 The values of m (degree) and n (order) as a function of the index
83 nm are given by the arrays indxm, indxn returned by getspecindx.
84
85 The associated legendre polynomials are normalized so that the
86 integral (pbar(n,m,theta)**2)*sin(theta) on the interval theta=0 to pi
87 is 1, where pbar(m,n,theta)=sqrt((2*n+1)*factorial(n-m)/(2*factorial(n+m)))*
88 sin(theta)**m/(2**n*factorial(n)) times the (n+m)th derivative of
89 (x**2-1)**n with respect to x=cos(theta).
90 theta = pi/2 - phi, where phi is latitude and theta is colatitude.
91 Therefore, cos(theta) = sin(phi) and sin(theta) = cos(phi).
92 Note that pbar(0,0,theta)=sqrt(2)/2, and pbar(1,0,theta)=.5*sqrt(6)*sin(lat).
93
94 The default grid type is regular (equally spaced latitude points).
95 Set gridtype='gaussian' when creating a class instance
96 for gaussian latitude points.
97
98 Quantities needed to compute spherical harmonics are precomputed and stored
99 when the class instance is created with legfunc='stored' (the default).
100 If legfunc='computed', they are recomputed on the fly on each method call.
101 The storage requirements for legfunc="stored" increase like nlat**2, while
102 those for legfunc='stored' increase like nlat**3. However, for
103 repeated method invocations on a single class instance, legfunc="stored"
104 will always be faster.
105
106 @contact: U{Jeff Whitaker<mailto:jeffrey.s.whitaker@noaa.gov>}
107
108 @version: 1.0.7
109
110 @license: Permission to use, copy, modify, and distribute this software and its
111 documentation for any purpose and without fee is hereby granted,
112 provided that the above copyright notice appear in all copies and that
113 both that copyright notice and this permission notice appear in
114 supporting documentation.
115 THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
116 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO
117 EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT OR
118 CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
119 USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
120 OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
121 PERFORMANCE OF THIS SOFTWARE.
122 """
123 import _spherepack, numpy, math, sys
124
125
126
127 _private_vars = ['nlon','nlat','gridtype','legfunc','rsphere']
128 __version__ = '1.0.9'
129
131 """
132 spherical harmonic transform class.
133
134 @ivar nlat: number of latitudes (set when class instance is created,
135 cannot be changed).
136
137 @ivar nlon: number of longitudes (set when class instance is created,
138 cannot be changed).
139
140 @ivar rsphere: The radius of the sphere in meters (set when class
141 instance is created, cannot be changed).
142
143 @ivar legfunc: 'stored' or 'computed'. If 'stored',
144 associated legendre functions are precomputed and stored when the
145 class instance is created. If 'computed', associated
146 legendre functions are computed on the fly when transforms are
147 requested. Set when class instance is created, cannot be changed.
148
149 @ivar gridtype: 'regular' (equally spaced in longitude and latitude)
150 or 'gaussian' (equally spaced in longitude, latitudes located at
151 roots of ordinary Legendre polynomial of degree nlat). Set when class
152 instance is created, cannot be changed.
153 """
154
156 """
157 prevent modification of read-only instance variables.
158 """
159 if key in self.__dict__ and key in _private_vars:
160 raise AttributeError('Attempt to rebind read-only instance variable '+key)
161 else:
162 self.__dict__[key] = val
163
165 """
166 prevent deletion of read-only instance variables.
167 """
168 if key in self.__dict__ and key in _private_vars:
169 raise AttributeError('Attempt to unbind read-only instance variable '+key)
170 else:
171 del self.__dict__[key]
172
173 - def __init__(self, nlon, nlat, rsphere=6.3712e6, gridtype='regular', legfunc='stored'):
174 """
175 create a Spharmt class instance.
176
177 @param nlon: Number of longitudes. The grid must be oriented from
178 east to west, with the first point at the Greenwich meridian
179 and the last point at 360-delta degrees east
180 (where delta = 360/nlon degrees). Must be >= 4. Transforms will
181 be faster when nlon is the product of small primes.
182
183 @param nlat: Number of latitudes. The grid must be oriented from north
184 to south. If nlat is odd the equator is included.
185 If nlat is even the equator will lie half way between points
186 points nlat/2 and (nlat/2)+1. Must be >=3.
187
188 @keyword rsphere: The radius of the sphere in meters.
189 Default 6371200 (the value for Earth).
190
191 @keyword legfunc: 'stored' (default) or 'computed'. If 'stored',
192 associated legendre functions are precomputed and stored when the
193 class instance is created. This uses O(nlat**3) memory, but
194 speeds up the spectral transforms. If 'computed', associated
195 legendre functions are computed on the fly when transforms are
196 requested. This uses O(nlat**2) memory, but slows down the spectral
197 transforms a bit.
198
199 @keyword gridtype: 'regular' (default) or 'gaussian'. Regular grids
200 will include the poles and equator if nlat is odd. Gaussian
201 grids never include the poles, but will include the equator if
202 nlat is odd.
203
204 """
205
206 if rsphere > 0.0:
207 self.rsphere= rsphere
208 else:
209 msg = 'Spharmt.__init__ illegal value of rsphere (%s) - must be postitive' % (rsphere)
210 raise ValueError(msg)
211 if nlon > 3:
212 self.nlon = nlon
213 else:
214 msg = 'Spharmt.__init__ illegal value of nlon (%s) - must be at least 4' % (nlon,)
215 raise ValueError(msg)
216 if nlat > 2:
217 self.nlat = nlat
218 else:
219 msg = 'Spharmt.__init__ illegal value of nlat (%s) - must be at least 3' % (nlat,)
220 raise ValueError(msg)
221 if gridtype != 'regular' and gridtype != 'gaussian':
222 msg = 'Spharmt.__init__ illegal value of gridtype (%s) - must be either "gaussian" or "regular"' % gridtype
223 raise ValueError(msg)
224 else:
225 self.gridtype = gridtype
226
227 if legfunc != 'computed' and legfunc != 'stored':
228 msg = 'Spharmt.__init__ illegal value of legfunc (%s) - must be either "computed" or "stored"' % legfunc
229 raise ValueError(msg)
230 else:
231 self.legfunc = legfunc
232
233 if nlon%2:
234 n1 = min(nlat, (nlon + 1)/2)
235 else:
236 n1 = min(nlat, (nlon + 2)/2)
237 if nlat%2:
238 n2 = (nlat + 1)/2
239 else:
240 n2 = nlat/2
241
242 if gridtype == 'regular':
243 if legfunc == 'stored':
244 lshaes = (n1*n2*(nlat + nlat - n1 + 1))/2 + nlon + 15
245 lwork = 5*nlat*n2 + 3*((n1 - 2)*(nlat + nlat - n1 -1))/2
246 wshaes, ierror = _spherepack.shaesi(nlat, nlon, lshaes, lwork, nlat+1)
247 if ierror != 0:
248 msg = 'In return from call to shaesi in Spharmt.__init__ ierror = %d' % ierror
249 raise ValueError(msg)
250 self.wshaes = wshaes
251 lshses = lshaes
252 wshses, ierror = _spherepack.shsesi(nlat, nlon, lshses, lwork, nlat+1)
253 if ierror != 0:
254 msg = 'In return from call to shsesi in Spharmt.__init__ ierror = %d' % ierror
255 raise ValueError(msg)
256 self.wshses = wshses
257 lvhaes = n1*n2*(nlat + nlat - n1 + 1) + nlon + 15
258 lwork = 3*(max(n1 -2,0)*(nlat + nlat - n1 - 1))/2 + 5*n2*nlat
259 wvhaes, ierror = _spherepack.vhaesi(nlat, nlon, lvhaes, lwork, 2*(nlat+1))
260 if ierror != 0:
261 msg = 'In return from call to vhaesi in Spharmt.__init__ ierror = %d' % ierror
262 raise ValueError(msg)
263 self.wvhaes = wvhaes
264 lwork = 3*(max(n1 - 2,0)*(nlat + nlat - n1 -1))/2 + 5*n2*nlat
265 lvhses = n1*n2*(nlat + nlat - n1 + 1) + nlon + 15
266 wvhses, ierror = _spherepack.vhsesi(nlat,nlon,lvhses,lwork,2*(nlat+1))
267 if ierror != 0:
268 msg = 'In return from call to vhsesi in Spharmt.__init__ ierror = %d' % ierror
269 raise ValueError(msg)
270 self.wvhses = wvhses
271 else:
272 lshaec = 2*nlat*n2+3*((n1-2)*(nlat+nlat-n1-1))/2+nlon+15
273 wshaec, ierror = _spherepack.shaeci(nlat, nlon, lshaec, 2*(nlat+1))
274 if ierror != 0:
275 msg = 'In return from call to shaeci in Spharmt.__init__ ierror = %d' % ierror
276 raise ValueError(msg)
277 self.wshaec = wshaec
278 lshsec = lshaec
279 wshsec, ierror = _spherepack.shseci(nlat, nlon, lshsec, 2*(nlat+1))
280 if ierror != 0:
281 msg = 'In return from call to shseci in Spharmt.__init__ ierror = %d' % ierror
282 raise ValueError(msg)
283 self.wshsec = wshsec
284 lvhaec = 4*nlat*n2+3*max(n1-2,0)*(2*nlat-n1-1)+nlon+15
285 wvhaec, ierror = _spherepack.vhaeci(nlat, nlon, lvhaec, 2*(nlat+1))
286 if ierror != 0:
287 msg = 'In return from call to vhaeci in Spharmt.__init__ ierror = %d' % ierror
288 raise ValueError(msg)
289 self.wvhaec = wvhaec
290 lvhsec = lvhaec
291 wvhsec, ierror = _spherepack.vhseci(nlat, nlon, lvhsec, 2*(nlat+1))
292 if ierror != 0:
293 msg = 'In return from call to vhseci in Spharmt.__init__ ierror = %d' % ierror
294 raise ValueError(msg)
295 self.wvhsec = wvhsec
296
297
298 elif gridtype == 'gaussian':
299 if legfunc == 'stored':
300 lshags = nlat*(3*(n1 + n2) - 2) + (n1 - 1)*(n2*(2*nlat - n1) - 3*n1)/2 + nlon + 15
301 lwork = 4*nlat*(nlat + 2) + 2
302 ldwork = nlat*(nlat + 4)
303 wshags, ierror = _spherepack.shagsi(nlat, nlon, lshags, lwork, ldwork)
304 if ierror != 0:
305 msg = 'In return from call to shagsi in Spharmt.__init__ ierror = %d' % ierror
306 raise ValueError(msg)
307 self.wshags = wshags
308 lshsgs = lshags
309 wshsgs, ierror = _spherepack.shsgsi(nlat, nlon, lshsgs, lwork, ldwork)
310 if ierror != 0:
311 msg = 'In return from call to shsgsi in Spharmt.__init__ ierror = %d' % ierror
312 raise ValueError(msg)
313 self.wshsgs = wshsgs
314 lvhags = (nlat +1)*(nlat + 1)*nlat/2 +nlon + 15
315 ldwork = (3*nlat*(nlat + 3) + 2)/2
316 wvhags, ierror = _spherepack.vhagsi(nlat, nlon, lvhags, ldwork)
317 if ierror != 0:
318 msg = 'In return from call to vhagsi in Spharmt.__init__ ierror = %d' % ierror
319 raise ValueError(msg)
320 self.wvhags = wvhags
321 lvhsgs = n1*n2*(nlat + nlat - n1 +1) + nlon + 15 + 2*nlat
322 ldwork = (3*nlat*(nlat + 3) + 2)/2
323 wvhsgs, ierror = _spherepack.vhsgsi(nlat, nlon, lvhsgs, ldwork)
324 if ierror != 0:
325 msg = 'In return from call to vhsgsi in Spharmt.__init__ ierror = %d' % ierror
326 raise ValueError(msg)
327 self.wvhsgs = wvhsgs
328 else:
329 lshagc = nlat*(2*n2+3*n1-2)+3*n1*(1-n1)/2+nlon+15
330 wshagc, ierror = _spherepack.shagci(nlat, nlon, lshagc, nlat*(nlat+4))
331 if ierror != 0:
332 msg = 'In return from call to shagci in Spharmt.__init__ ierror = %d' % ierror
333 raise ValueError(msg)
334 self.wshagc = wshagc
335 lshsgc = lshagc
336 wshsgc, ierror = _spherepack.shsgci(nlat, nlon, lshsgc, nlat*(nlat+4))
337 if ierror != 0:
338 msg = 'In return from call to shsgci in Spharmt.__init__ ierror = %d' % ierror
339 raise ValueError(msg)
340 self.wshsgc = wshsgc
341 lvhagc = 4*nlat*n2+3*max(n1-2,0)*(2*nlat-n1-1)+nlon+n2+15
342 ldwork = 2*nlat*(nlat+1)+1
343 wvhagc, ierror = _spherepack.vhagci(nlat, nlon, lvhagc, ldwork)
344 if ierror != 0:
345 msg = 'In return from call to vhagci in Spharmt.__init__ ierror = %d' % ierror
346 raise ValueError(msg)
347 self.wvhagc = wvhagc
348 lvhsgc = 4*nlat*n2+3*max(n1-2,0)*(2*nlat-n1-1)+nlon+15
349 wvhsgc, ierror = _spherepack.vhsgci(nlat, nlon, lvhsgc, ldwork)
350 if ierror != 0:
351 msg = 'In return from call to vhsgci in Spharmt.__init__ ierror = %d' % ierror
352 raise ValueError(msg)
353 self.wvhsgc = wvhsgc
354
356 """
357 grid to spectral transform (spherical harmonic analysis).
358
359 @param datagrid: rank 2 or 3 numpy float32 array with shape (nlat,nlon) or
360 (nlat,nlon,nt), where nt is the number of grids to be transformed. If
361 datagrid is rank 2, nt is assumed to be 1.
362
363 @keyword ntrunc: optional spectral truncation limit.
364 (default self.nlat-1)
365
366 @return: C{B{dataspec}} - rank 1 or 2 numpy complex array with shape
367 (ntrunc+1)*(ntrunc+2)/2 or ((ntrunc+1)*(ntrunc+2)/2,nt) containing
368 complex spherical harmonic coefficients resulting from the spherical
369 harmonic analysis of datagrid.
370 """
371
372
373
374
375 idim = datagrid.ndim
376
377 if len(datagrid.shape) > 3:
378 msg = 'grdtospec needs a rank two or three array, got %d' % (len(datagrid.shape),)
379 raise ValueError(msg)
380
381 if datagrid.shape[0] != self.nlat or datagrid.shape[1] != self.nlon:
382 msg = 'grdtospec needs an array of size %d by %d, got %d by %d' % (self.nlat, self.nlon, datagrid.shape[0], datagrid.shape[1],)
383 raise ValueError(msg)
384
385
386
387 if ntrunc is None:
388 ntrunc = self.nlat-1
389
390 if ntrunc < 0 or ntrunc+1 > datagrid.shape[0]:
391 msg = 'ntrunc must be between 0 and %d' % (datagrid.shape[0]-1,)
392 raise ValueError(msg)
393
394
395 nlat = self.nlat
396 nlon = self.nlon
397 if nlat%2:
398 n2 = (nlat + 1)/2
399 else:
400 n2 = nlat/2
401
402 if len(datagrid.shape) == 2:
403 nt = 1
404 datagrid = numpy.reshape(datagrid, (nlat,nlon,1))
405 else:
406 nt = datagrid.shape[2]
407
408 if self.gridtype == 'regular':
409
410
411
412 if self.legfunc == 'stored':
413 lwork = (nt+1)*nlat*nlon
414 a,b,ierror = _spherepack.shaes(datagrid,self.wshaes,lwork)
415 if ierror != 0:
416 msg = 'In return from call to shaes in Spharmt.grdtospec ierror = %d' % ierror
417 raise ValueError(msg)
418 else:
419 lwork = nlat*(nt*nlon+max(3*n2,nlon))
420 a,b,ierror = _spherepack.shaec(datagrid,self.wshaec,lwork)
421 if ierror != 0:
422 msg = 'In return from call to shaec in Spharmt.grdtospec ierror = %d' % ierror
423
424 raise ValueError(msg)
425
426
427
428 elif self.gridtype == 'gaussian':
429
430
431
432 if self.legfunc == 'stored':
433 lwork = nlat*nlon*(nt+1)
434 a,b,ierror = _spherepack.shags(datagrid,self.wshags,lwork)
435 if ierror != 0:
436 msg = 'In return from call to shags in Spharmt.grdtospec ierror = %d' % ierror
437 raise ValueError(msg)
438 else:
439 lwork = nlat*(nlon*nt+max(3*n2,nlon))
440 a,b,ierror = _spherepack.shagc(datagrid,self.wshagc,lwork)
441 if ierror != 0:
442 msg = 'In return from call to shagc in Spharmt.grdtospec ierror = %d' % ierror
443 raise ValueError(msg)
444
445
446
447 dataspec = _spherepack.twodtooned(a,b,ntrunc)
448
449 if idim == 2:
450 return numpy.squeeze(dataspec)
451 else:
452 return dataspec
453
455
456 """
457 spectral to grid transform (spherical harmonic synthesis).
458
459 @param dataspec: rank 1 or 2 numpy complex array with shape
460 (ntrunc+1)*(ntrunc+2)/2 or ((ntrunc+1)*(ntrunc+2)/2,nt) containing
461 complex spherical harmonic coefficients (where ntrunc is the
462 triangular truncation limit and nt is the number of spectral arrays
463 to be transformed). If dataspec is rank 1, nt is assumed to be 1.
464
465 @return: C{B{datagrid}} - rank 2 or 3 numpy float32 array with shape
466 (nlat,nlon) or (nlat,nlon,nt) containing the gridded data resulting from
467 the spherical harmonic synthesis of dataspec.
468 """
469
470
471
472 idim = dataspec.ndim
473
474 if len(dataspec.shape) > 2:
475 msg = 'spectogrd needs a rank one or two array, got %d' % (len(dataspec.shape),)
476 raise ValueError(msg)
477
478 nlat = self.nlat
479 nlon = self.nlon
480 if nlat%2:
481 n2 = (nlat + 1)/2
482 else:
483 n2 = nlat/2
484
485 if len(dataspec.shape) == 1:
486 nt = 1
487 dataspec = numpy.reshape(dataspec, (dataspec.shape[0],1))
488 else:
489 nt = dataspec.shape[1]
490
491 ntrunc = int(-1.5 + 0.5*math.sqrt(9.-8.*(1.-dataspec.shape[0])))
492 if ntrunc > nlat-1:
493 msg = 'ntrunc too large - can be max of %d, got %d' % (nlat-1,ntrunc)
494 raise ValueError(msg)
495
496 a, b = _spherepack.onedtotwod(dataspec,nlat)
497
498
499
500 if self.gridtype == 'regular':
501
502
503
504 if self.legfunc == 'stored':
505 lwork = (nt+1)*nlat*nlon
506 datagrid, ierror = _spherepack.shses(nlon,a,b,self.wshses,lwork)
507 if ierror != 0:
508 msg = 'In return from call to shses in Spharmt.spectogrd ierror = %d' % ierror
509 raise ValueError(msg)
510 else:
511 lwork = nlat*(nt*nlon+max(3*n2,nlon))
512 datagrid, ierror = _spherepack.shsec(nlon,a,b,self.wshsec,lwork)
513 if ierror != 0:
514 msg = 'In return from call to shsec in Spharmt.spectogrd ierror = %d' % ierror
515 raise ValueError(msg)
516
517
518
519 elif self.gridtype == 'gaussian':
520
521
522
523 if self.legfunc == 'stored':
524 lwork = nlat*nlon*(nt+1)
525 datagrid, ierror = _spherepack.shsgs(nlon,a,b,self.wshsgs,lwork)
526 if ierror != 0:
527 msg = 'In return from call to shsgs in Spharmt.spectogrd ierror = %d' % ierror
528 raise ValueError(msg)
529 else:
530 lwork = nlat*(nlon*nt+max(3*n2,nlon))
531 datagrid, ierror = _spherepack.shsgc(nlon,a,b,self.wshsgc,lwork)
532 if ierror != 0:
533 msg = 'In return from call to shsgc in Spharmt.spectogrd ierror = %d' % ierror
534 raise ValueError(msg)
535
536 if idim == 1:
537 return numpy.squeeze(datagrid)
538 else:
539 return datagrid
540
542
543 """
544 compute spectral coefficients of vorticity and divergence given vector wind.
545
546 @param ugrid: rank 2 or 3 numpy float32 array containing grid of zonal
547 winds. Must have shape (nlat,nlon) or (nlat,nlon,nt), where nt is the number
548 of grids to be transformed. If ugrid is rank 2, nt is assumed to be 1.
549
550 @param vgrid: rank 2 or 3 numpy float32 array containing grid of meridional
551 winds. Must have shape (nlat,nlon) or (nlat,nlon,nt), where nt is the number
552 of grids to be transformed. Both ugrid and vgrid must have the same shape.
553
554 @keyword ntrunc: optional spectral truncation limit.
555 (default self.nlat-1)
556
557 @return: C{B{vrtspec, divspec}} - rank 1 or 2 numpy complex arrays
558 of vorticity and divergence spherical harmonic coefficients with shape
559 shape (ntrunc+1)*(ntrunc+2)/2 or ((ntrunc+1)*(ntrunc+2)/2,nt).
560 """
561
562
563
564 idim = ugrid.ndim
565
566 shapeu = ugrid.shape
567 shapev = vgrid.shape
568
569 if ntrunc is None:
570 ntrunc = self.nlat-1
571
572 if shapeu != shapev:
573 msg = 'getvrtdivspec input arrays must be same shape!'
574 raise ValueError(msg)
575
576
577 if len(shapeu) !=2 and len(shapeu) !=3:
578 msg = 'getvrtdivspec needs rank two or three arrays!'
579 raise ValueError(msg)
580
581 if shapeu[0] != self.nlat or shapeu[1] != self.nlon:
582 msg = 'getvrtdivspec needs input arrays whose first two dimensions are si%d and %d, got %d and %d' % (self.nlat, self.nlon, ugrid.shape[0], ugrid.shape[1],)
583 raise ValueError(msg)
584
585
586
587
588 if ntrunc < 0 or ntrunc+1 > shapeu[0]:
589 msg = 'ntrunc must be between 0 and %d' % (ugrid.shape[0]-1,)
590 raise ValueError(msg)
591
592 nlat = self.nlat
593 nlon = self.nlon
594 if nlat%2:
595 n2 = (nlat + 1)/2
596 else:
597 n2 = nlat/2
598 rsphere= self.rsphere
599
600
601
602
603 if len(shapeu) == 2:
604 nt = 1
605 w = numpy.reshape(ugrid, (nlat,nlon,1))
606 v = -numpy.reshape(vgrid, (nlat,nlon,1))
607 else:
608 nt = shapeu[2]
609 w = ugrid
610 v = -vgrid
611
612
613
614 if self.gridtype == 'regular':
615
616
617
618 if self.legfunc == 'stored':
619 lwork = (2*nt+1)*nlat*nlon
620 br,bi,cr,ci,ierror = _spherepack.vhaes(v,w,self.wvhaes,lwork)
621 if ierror != 0:
622 msg = 'In return from call to vhaes in Spharmt.getvrtdivspec ierror = %d' % ierror
623 raise ValueError(msg)
624 else:
625 lwork = nlat*(2*nt*nlon+max(6*n2,nlon))
626 br,bi,cr,ci,ierror = _spherepack.vhaec(v,w,self.wvhaec,lwork)
627 if ierror != 0:
628 msg = 'In return from call to vhaec in Spharmt.getvrtdivspec ierror = %d' % ierror
629 raise ValueError(msg)
630
631
632
633 elif self.gridtype == 'gaussian':
634
635
636
637 if self.legfunc == 'stored':
638 lwork = (2*nt+1)*nlat*nlon
639 br,bi,cr,ci,ierror = _spherepack.vhags(v,w,self.wvhags,lwork)
640 if ierror != 0:
641 msg = 'In return from call to vhags in Spharmt.getvrtdivspec ierror = %d' % ierror
642 raise ValueError(msg)
643 else:
644 lwork = 2*nlat*(2*nlon*nt+3*n2)
645 br,bi,cr,ci,ierror = _spherepack.vhagc(v,w,self.wvhagc,lwork)
646 if ierror != 0:
647 msg = 'In return from call to vhagc in Spharmt.getvrtdivspec ierror = %d' % ierror
648 raise ValueError(msg)
649
650
651
652
653 vrtspec, divspec = _spherepack.twodtooned_vrtdiv(br,bi,cr,ci,ntrunc,rsphere)
654
655 if idim == 2:
656 return numpy.squeeze(vrtspec), numpy.squeeze(divspec)
657 else:
658 return vrtspec, divspec
659
660 - def getuv(self, vrtspec, divspec):
661
662 """
663 compute vector wind on grid given complex spectral coefficients
664 of vorticity and divergence.
665
666 @param vrtspec: rank 1 or 2 numpy complex array of vorticity spectral
667 coefficients, with shape (ntrunc+1)*(ntrunc+2)/2 or
668 ((ntrunc+1)*(ntrunc+2)/2,nt) (where ntrunc is the triangular truncation
669 and nt is the number of spectral arrays to be transformed).
670 If vrtspec is rank 1, nt is assumed to be 1.
671
672 @param divspec: rank 1 or 2 numpy complex array of divergence spectral
673 coefficients, with shape (ntrunc+1)*(ntrunc+2)/2 or
674 ((ntrunc+1)*(ntrunc+2)/2,nt) (where ntrunc is the triangular truncation
675 and nt is the number of spectral arrays to be transformed).
676 Both vrtspec and divspec must have the same shape.
677
678 @return: C{B{ugrid, vgrid}} - rank 2 or 3 numpy float32 arrays containing
679 gridded zonal and meridional winds. Shapes are either (nlat,nlon) or
680 (nlat,nlon,nt).
681 """
682
683 idim = vrtspec.ndim
684
685 shapevrt = vrtspec.shape
686 shapediv = divspec.shape
687
688
689
690 if shapevrt != shapediv:
691 msg = 'vrtspec, divspec must be same size in getuv!'
692 raise ValueError(msg)
693
694 if len(shapevrt) !=1 and len(shapevrt) !=2:
695 msg = 'getuv needs rank one or two input arrays!'
696 raise ValueError(msg)
697
698
699
700
701 ntrunc = int(-1.5 + 0.5*math.sqrt(9.-8.*(1.-vrtspec.shape[0])))
702
703 nlat = self.nlat
704 nlon = self.nlon
705 if nlat%2:
706 n2 = (nlat + 1)/2
707 else:
708 n2 = nlat/2
709 rsphere= self.rsphere
710
711 if len(vrtspec.shape) == 1:
712 nt = 1
713 vrtspec = numpy.reshape(vrtspec, (vrtspec.shape[0],1))
714 divspec = numpy.reshape(divspec, (divspec.shape[0],1))
715 else:
716 nt = vrtspec.shape[1]
717
718
719
720 br,bi,cr,ci = _spherepack.onedtotwod_vrtdiv(vrtspec,divspec,nlat,rsphere)
721
722
723
724 if self.gridtype == 'regular':
725
726
727
728 if self.legfunc == 'stored':
729 lwork = (2*nt+1)*nlat*nlon
730 v, w, ierror = _spherepack.vhses(nlon,br,bi,cr,ci,self.wvhses,lwork)
731 if ierror != 0:
732 msg = 'In return from call to vhses in Spharmt.getuv ierror = %d' % ierror
733 raise ValueError(msg)
734 else:
735 lwork = nlat*(2*nt*nlon+max(6*n2,nlon))
736 v, w, ierror = _spherepack.vhsec(nlon,br,bi,cr,ci,self.wvhsec,lwork)
737 if ierror != 0:
738 msg = 'In return from call to vhsec in Spharmt.getuv ierror = %d' % ierror
739 raise ValueError(msg)
740
741
742
743 elif self.gridtype == 'gaussian':
744
745
746
747 if self.legfunc == 'stored':
748 lwork = (2*nt+1)*nlat*nlon
749 v, w, ierror = _spherepack.vhsgs(nlon,br,bi,cr,ci,self.wvhsgs,lwork)
750 if ierror != 0:
751 msg = 'In return from call to vhsgs in Spharmt.getuv ierror = %d' % ierror
752 raise ValueError(msg)
753 else:
754 lwork = nlat*(2*nt*nlon+max(6*n2,nlon))
755 v, w, ierror = _spherepack.vhsgc(nlon,br,bi,cr,ci,self.wvhsgc,lwork)
756 if ierror != 0:
757 msg = 'In return from call to vhsgc in Spharmt.getuv ierror = %d' % ierror
758 raise ValueError(msg)
759
760
761
762 if idim == 1:
763 return numpy.reshape(w, (nlat,nlon)), -numpy.reshape(v, (nlat,nlon))
764 else:
765 return w,-v
766
767 - def getpsichi(self, ugrid, vgrid, ntrunc=None):
768
769 """
770 compute streamfunction and velocity potential on grid given vector wind.
771
772 @param ugrid: rank 2 or 3 numpy float32 array containing grid of zonal
773 winds. Must have shape (nlat,nlon) or (nlat,nlon,nt), where nt is the number
774 of grids to be transformed. If ugrid is rank 2, nt is assumed to be 1.
775
776 @param vgrid: rank 2 or 3 numpy float32 array containing grid of meridional
777 winds. Must have shape (nlat,nlon) or (nlat,nlon,nt), where nt is the number
778 of grids to be transformed. Both ugrid and vgrid must have the same shape.
779
780 @keyword ntrunc: optional spectral truncation limit.
781 (default self.nlat-1)
782
783 @return: C{B{psigrid, chigrid}} - rank 2 or 3 numpy float32 arrays
784 of gridded streamfunction and velocity potential. Shapes are either
785 (nlat,nlon) or (nlat,nlon,nt).
786 """
787
788
789
790 idim = ugrid.ndim
791
792 shapeu = ugrid.shape
793 shapev = vgrid.shape
794
795 if ntrunc is None:
796 ntrunc = self.nlat-1
797
798 if shapeu != shapev:
799 msg = 'getvrtdivspec input arrays must be same shape!'
800 raise ValueError(msg)
801
802
803 if len(shapeu) !=2 and len(shapeu) !=3:
804 msg = 'getvrtdivspec needs rank two or three arrays!'
805 raise ValueError(msg)
806
807 if shapeu[0] != self.nlat or shapeu[1] != self.nlon:
808 msg = 'getpsichi needs input arrays whose first two dimensions are si%d and %d, got %d and %d' % (self.nlat, self.nlon, ugrid.shape[0], ugrid.shape[1],)
809 raise ValueError(msg)
810
811
812
813 if ntrunc < 0 or ntrunc+1 > ugrid.shape[0]:
814 msg = 'ntrunc must be between 0 and %d' % (ugrid.shape[0]-1,)
815 raise ValueError(msg)
816
817
818
819 vrtspec, divspec = self.getvrtdivspec(ugrid, vgrid, ntrunc)
820
821
822
823 if len(vrtspec.shape) == 1:
824 nt = 1
825 vrtspec = numpy.reshape(vrtspec, ((ntrunc+1)*(ntrunc+2)/2,1))
826 divspec = numpy.reshape(divspec, ((ntrunc+1)*(ntrunc+2)/2,1))
827 else:
828 nt = vrtspec.shape[1]
829
830
831
832 psispec = _spherepack.invlap(vrtspec, self.rsphere)
833 chispec = _spherepack.invlap(divspec, self.rsphere)
834
835
836
837 psigrid = self.spectogrd(psispec)
838 chigrid = self.spectogrd(chispec)
839
840 if idim == 2:
841 return numpy.squeeze(psigrid), numpy.squeeze(chigrid)
842 else:
843 return psigrid, chigrid
844
846
847 """
848 compute vector gradient on grid given complex spectral coefficients.
849
850 @param chispec: rank 1 or 2 numpy complex array with shape
851 (ntrunc+1)*(ntrunc+2)/2 or ((ntrunc+1)*(ntrunc+2)/2,nt) containing
852 complex spherical harmonic coefficients (where ntrunc is the
853 triangular truncation limit and nt is the number of spectral arrays
854 to be transformed). If chispec is rank 1, nt is assumed to be 1.
855
856 @return: C{B{uchi, vchi}} - rank 2 or 3 numpy float32 arrays containing
857 gridded zonal and meridional components of the vector gradient.
858 Shapes are either (nlat,nlon) or (nlat,nlon,nt).
859 """
860
861
862
863 idim = chispec.ndim
864
865 if len(chispec.shape) !=1 and len(chispec.shape) !=2:
866 msg = 'getgrad needs rank one or two arrays!'
867 raise ValueError(msg)
868
869
870
871
872 ntrunc = int(-1.5 + 0.5*math.sqrt(9.-8.*(1.-chispec.shape[0])))
873
874
875
876 if len(chispec.shape) == 1:
877 nt = 1
878 chispec = numpy.reshape(chispec, ((ntrunc+1)*(ntrunc+2)/2,1))
879 else:
880 nt = chispec.shape[1]
881
882
883
884 divspec = _spherepack.lap(chispec,self.rsphere)
885
886
887
888 uchi, vchi = self.getuv(numpy.zeros(chispec.shape, chispec.dtype), divspec)
889
890 if idim == 1:
891 return numpy.squeeze(uchi), numpy.squeeze(vchi)
892 else:
893 return uchi, vchi
894
896
897 """
898 isotropic spectral smoothing on a sphere.
899
900 @param datagrid: rank 2 or 3 numpy float32 array with shape (nlat,nlon) or
901 (nlat,nlon,nt), where nt is the number of grids to be smoothed. If
902 datagrid is rank 2, nt is assumed to be 1.
903
904 @param smooth: rank 1 array of length nlat containing smoothing factors
905 as a function of total wavenumber.
906
907 @return: C{B{datagrid}} - rank 2 or 3 numpy float32 array with shape
908 (nlat,nlon) or (nlat,nlon,nt) containing the smoothed grids.
909 """
910
911
912
913
914 if len(datagrid.shape) > 3:
915 msg = 'specsmooth needs a rank two or three array, got %d' % (len(datagrid.shape),)
916 raise ValueError(msg)
917
918 if datagrid.shape[0] != self.nlat or datagrid.shape[1] != self.nlon:
919 msg = 'specsmooth needs an array of size %d by %d, got %d by %d' % (self.nlat, self.nlon, datagrid.shape[0], datagrid.shape[1],)
920 raise ValueError(msg)
921
922
923
924
925 if len(smooth.shape) !=1 or smooth.shape[0] != datagrid.shape[0]:
926 msg = 'smooth must be rank 1 and same size as datagrid.shape[0] in specsmooth!'
927 raise ValueError(msg)
928
929
930
931 nlat = self.nlat
932 dataspec = self.grdtospec(datagrid, nlat-1)
933
934
935
936 dataspec = _spherepack.multsmoothfact(dataspec, smooth)
937
938
939
940 datagrid = self.spectogrd(dataspec)
941
942 return datagrid
943
944 -def regrid(grdin, grdout, datagrid, ntrunc=None, smooth=None):
945 """
946 regrid data using spectral interpolation, while performing
947 optional spectral smoothing and/or truncation.
948
949 @param grdin: Spharmt class instance describing input grid.
950
951 @param grdout: Spharmt class instance describing output grid.
952
953 @param datagrid: data on input grid (grdin.nlat x grdin.nlon). If
954 datagrid is rank 3, last dimension is the number of grids to interpolate.
955
956 @keyword ntrunc: optional spectral truncation limit for datagrid
957 (default min(grdin.nlat-1,grdout.nlat-1)).
958
959 @keyword smooth: rank 1 array of length grdout.nlat containing smoothing
960 factors as a function of total wavenumber (default is no smoothing).
961
962 @return: C{B{datagrid}} - interpolated (and optionally smoothed) array(s)
963 on grdout.nlon x grdout.nlat grid.
964 """
965
966
967
968
969 if len(datagrid.shape) > 3:
970 msg = 'regrid needs a rank two or three array, got %d' % (len(datagrid.shape),)
971 raise ValueError(msg)
972
973 if datagrid.shape[0] != grdin.nlat or datagrid.shape[1] != grdin.nlon:
974 msg = 'grdtospec needs an array of size %d by %d, got %d by %d' % (grdin.nlat, grdin.nlon, datagrid.shape[0], datagrid.shape[1],)
975 raise ValueError(msg)
976
977 if smooth is not None and (len(smooth.shape) !=1 or smooth.shape[0] != grdout.nlat):
978 msg = 'smooth must be rank 1 size grdout.nlat in regrid!'
979 raise ValueError(msg)
980
981 if ntrunc is None:
982 ntrunc = min(grdout.nlat-1,grdin.nlat-1)
983
984 dataspec = grdin.grdtospec(datagrid,ntrunc)
985
986 if smooth is not None:
987 dataspec = _spherepack.multsmoothfact(dataspec, smooth)
988
989 return grdout.spectogrd(dataspec)
990
992
993 """
994 compute the gaussian latitudes (in degrees) and quadrature weights.
995
996 @param nlat: number of gaussian latitudes desired.
997
998 @return: C{B{lats, wts}} - rank 1 numpy float64 arrays containing
999 gaussian latitudes (in degrees north) and gaussian quadrature weights.
1000 """
1001
1002
1003
1004 colats, wts, ierror = _spherepack.gaqd(nlat)
1005
1006 if ierror != 0:
1007 msg = 'In return from call to gaqd ierror = %d' % ierror
1008 raise ValueError(msg)
1009
1010
1011
1012 lats = 90.0 - colats*180.0/math.pi
1013
1014 return lats, wts
1015
1017
1018 """
1019 compute indices of zonal wavenumber (indxm) and degree (indxn)
1020 for complex spherical harmonic coefficients.
1021
1022 @param ntrunc: spherical harmonic triangular truncation limit.
1023
1024 @return: C{B{indxm, indxn}} - rank 1 numpy Int32 arrays
1025 containing zonal wavenumber (indxm) and degree (indxn) of
1026 spherical harmonic coefficients.
1027 """
1028
1029 indexn = numpy.indices((ntrunc+1,ntrunc+1))[1,:,:]
1030 indexm = numpy.indices((ntrunc+1,ntrunc+1))[0,:,:]
1031 indices = numpy.nonzero(numpy.greater(indexn, indexm-1).flatten())
1032 indxn = numpy.take(indexn.flatten(),indices)
1033 indxm = numpy.take(indexm.flatten(),indices)
1034
1035 return numpy.squeeze(indxm), numpy.squeeze(indxn)
1036
1038 """
1039 computes the lat/lon values of the points on the surface of the sphere
1040 corresponding to a twenty-sided (icosahedral) geodesic.
1041
1042 @param m: the number of points on the edge of a single geodesic triangle.
1043 There are 10*(m-1)**2+2 total geodesic points, including the poles.
1044
1045 @return: C{B{lats, lons}} - rank 1 numpy float32 arrays containing
1046 the latitudes and longitudes of the geodesic points (in degrees). These
1047 points are nearly evenly distributed on the surface of the sphere.
1048 """
1049 x,y,z = _spherepack.ihgeod(m)
1050
1051 rad2dg = 180./math.pi
1052 r1 = x*x+y*y
1053 r = numpy.sqrt(r1+z*z)
1054 r1 = numpy.sqrt(r1)
1055 xtmp = numpy.where(numpy.logical_or(x,y),x,numpy.ones(x.shape,numpy.float32))
1056 ztmp = numpy.where(numpy.logical_or(r1,z),z,numpy.ones(z.shape,numpy.float32))
1057 lons = rad2dg*numpy.arctan2(y,xtmp)+180.
1058 lats = rad2dg*numpy.arctan2(r1,ztmp)-90.
1059 lat = numpy.zeros(10*(m-1)**2+2,numpy.float32)
1060 lon = numpy.zeros(10*(m-1)**2+2,numpy.float32)
1061
1062 lat[0] = 90; lat[1] = -90.
1063 lon[0] = 0.; lon[1] = 0.
1064 lat[2:] = lats[0:2*(m-1),0:m-1,:].flatten()
1065 lon[2:] = lons[0:2*(m-1),0:m-1,:].flatten()
1066 return lat,lon
1067
1069 """
1070 calculate associated legendre functions for triangular truncation T(ntrunc),
1071 at a given latitude.
1072
1073 @param lat: the latitude (in degrees) to compute the associate legendre
1074 functions.
1075
1076 @param ntrunc: the triangular truncation limit.
1077
1078 @return: C{B{pnm}} - rank 1 numpy float32 array containing the
1079 (C{B{ntrunc}}+1)*(C{B{ntrunc}}+2)/2 associated legendre functions at
1080 latitude C{B{lat}}.
1081 """
1082 return _spherepack.getlegfunc(lat,ntrunc)
1083
1085 """
1086 spectral interpolation given spherical harmonic coefficients.
1087
1088 @param lon: longitude (in degrees) of point on a sphere to interpolate to.
1089
1090 @param dataspec: spectral coefficients of function to interpolate.
1091
1092 @param legfuncs: associated legendre functions with same triangular
1093 truncation as C{B{dataspec}} (computed using L{legendre}), computed
1094 at latitude of interpolation point.
1095
1096 @return: C{B{ob}} - interpolated value.
1097 """
1098 ntrunc1 = int(-1.5 + 0.5*math.sqrt(9.-8.*(1.-dataspec.shape[0])))
1099 ntrunc2 = int(-1.5 + 0.5*math.sqrt(9.-8.*(1.-legfuncs.shape[0])))
1100 if ntrunc1 != ntrunc2:
1101 raise ValueError('first dimensions of dataspec and legfuncs in Spharmt.specintrp imply inconsistent spectral truncations - they must be the same!')
1102 return _spherepack.specintrp((math.pi/180.)*lon,ntrunc1,dataspec,legfuncs)
1103