VAPoR  0.1
vapor_utils.py
Go to the documentation of this file.
1 ''' The vapor_utils module contains:
2  mag3d - calculate magnitude of 3-vector
3  mag2d - calculate magnitude of 2-vector
4  deriv_findiff - calculate derivative using 6th order finite differences
5  curl_findiff - calculate curl using 6th order finite differences
6  div_findiff - calculate divergence using 6 order finite differences
7  grad_findiff - calculate gradient using 6th order finite differences.
8  deriv_var_findiff - calculate a derivative of one 3D
9 variable with respect to another variable.
10  interp3d - interpolate a 3D variable to a vertical level surface of another variable.
11  vector_rotate - rotate and scale vector field for lat-lon grid.
12 '''
13 
14 def mag3d(a1,a2,a3):
15  '''Calculate the magnitude of a 3-vector.
16  Calling sequence: MAG = mag3d(A,B,C)
17  Where: A, B, and C are float32 numpy arrays.
18  Result MAG is a float32 numpy array containing the square root
19  of the sum of the squares of A, B, and C.'''
20  from numpy import sqrt
21  return sqrt(a1*a1 + a2*a2 + a3*a3)
22 
23 def mag2d(a1,a2):
24  '''Calculate the magnitude of a 2-vector.
25  Calling sequence: MAG = mag2d(A,B)
26  Where: A, and B are float32 numpy arrays.
27  Result MAG is a float32 numpy array containing the square root
28  of the sum of the squares of A and B.'''
29  from numpy import sqrt
30  return sqrt(a1*a1 + a2*a2)
31 
32 def deriv_findiff2(a,dir,dx):
33  ''' Function that calculates first-order derivatives
34  using 2nd order finite differences in regular Cartesian grids.'''
35 
36  import numpy
37  s = numpy.shape(a) #size of the input array
38  aprime = numpy.array(a) #output has the same size than input
39 
40  #
41  # derivative for user dir=1, in python this is third coordinate
42  #
43  if dir == 1:
44 
45  #forward differences near the first boundary
46  for i in range(1):
47  aprime[:,:,i] = (-a[:,:,i]+a[:,:,i+1]) / (dx)
48 
49  #centered differences
50  for i in range(1,s[2]-1):
51  aprime[:,:,i] = (-a[:,:,i-1]+a[:,:,i+1])/(2*dx)
52 
53  #backward differences near the second boundary
54  for i in range(s[2]-1,s[2]):
55  aprime[:,:,i] = (a[:,:,i-1]-a[:,:,i]) /(dx)
56 
57  #
58  # derivative for dir=2
59  #
60  if dir == 2:
61  #forward differences near the first boundary
62  for i in range(1):
63  aprime[:,i,:] = (-a[:,i,:]+a[:,i+1,:]) / (dx)
64 
65  #centered differences
66  for i in range(1,s[1]-1):
67  aprime[:,i,:] = (-a[:,i-1,:]+a[:,i+1,:])/(2*dx)
68 
69  #backward differences near the second boundary
70  for i in range(s[1]-1,s[1]):
71  aprime[:,i,:] = (a[:,i-1,:]-a[:,i,:]) /(dx)
72 
73  #
74  # derivative for user dir=3
75  #
76  if dir == 3:
77  #forward differences near the first boundary
78  for i in range(1):
79  aprime[i,:,:] = (-a[i,:,:]+a[i+1,:,:]) / (dx)
80 
81  #centered differences
82  for i in range(1,s[0]-1):
83  aprime[i,:,:] = (-a[i-1,:,:]+a[i+1,:,:])/(2*dx)
84 
85  #backward differences near the second boundary
86  for i in range(s[0]-1,s[0]):
87  aprime[i,:,:] = (a[i-1,:,:]-a[i,:,:]) /(dx)
88 
89  return aprime
90 
91 def deriv_findiff4(a,dir,dx):
92  ''' Function that calculates first-order derivatives
93  using 4th order finite differences in regular Cartesian grids.'''
94 
95  import numpy
96  s = numpy.shape(a) #size of the input array
97  aprime = numpy.array(a) #output has the same size than input
98 
99  #
100  # derivative for user dir=1, in python this is third coordinate
101  #
102  if dir == 1:
103 
104  #forward differences near the first boundary
105  for i in range(2):
106  aprime[:,:,i] = (-3*a[:,:,i]+4*a[:,:,i+1]-a[:,:,i+2]) / (2*dx)
107 
108  #centered differences
109  for i in range(2,s[2]-2):
110  aprime[:,:,i] = (a[:,:,i-2]-8*a[:,:,i-1]+8*a[:,:,i+1]-a[:,:,i+2])/(12*dx)
111 
112  #backward differences near the second boundary
113  for i in range(s[2]-2,s[2]):
114  aprime[:,:,i] = (a[:,:,i-2]-4*a[:,:,i-1]+3*a[:,:,i]) /(2*dx)
115 
116  #
117  # derivative for dir=2
118  #
119  if dir == 2:
120  #forward differences near the first boundary
121  for i in range(2):
122  aprime[:,i,:] = (-3*a[:,i,:]+4*a[:,i+1,:]-a[:,i+2,:]) / (2*dx)
123 
124  #centered differences
125  for i in range(2,s[1]-2):
126  aprime[:,i,:] = (a[:,i-2,:]-8*a[:,i-1,:]+8*a[:,i+1,:]-a[:,i+2,:])/(12*dx)
127 
128  #backward differences near the second boundary
129  for i in range(s[1]-2,s[1]):
130  aprime[:,i,:] = (a[:,i-2,:]-4*a[:,i-1,:]+3*a[:,i,:]) /(2*dx)
131 
132  #
133  # derivative for user dir=3
134  #
135  if dir == 3:
136  #forward differences near the first boundary
137  for i in range(2):
138  aprime[i,:,:] = (-3*a[i,:,:]+4*a[i+1,:,:]-a[i+2,:,:]) / (2*dx)
139 
140  #centered differences
141  for i in range(2,s[0]-2):
142  aprime[i,:,:] = (a[i-2,:,:]-8*a[i-1,:,:]+8*a[i+1,:,:]-a[i+2,:,:])/(12*dx)
143 
144  #backward differences near the second boundary
145  for i in range(s[0]-2,s[0]):
146  aprime[i,:,:] = (a[i-2,:,:]-4*a[i-1,:,:]+3*a[i,:,:]) /(2*dx)
147 
148 
149  return aprime
150 
151 def deriv_findiff(a,dir,dx,order=6):
152  ''' Function that calculates first-order derivatives
153  using sixth order finite differences in regular Cartesian grids.
154  Calling sequence: deriv = deriv_findiff(ary,dir,delta, order=6)
155  ary is a 3-dimensional float32 numpy array
156  dir is 1, 2, or 3 (for X, Y, or Z directions in user coordinates)
157  delta is the grid spacing in the direction dir.
158  order is the accuracy order (one of 6,4,2)
159  Returned array 'deriv' is the derivative of ary in the
160  specified direction dir. '''
161  import numpy
162 
163  if order == 4:
164  return deriv_findiff4(a,dir,dx)
165 
166  if order == 2:
167  return deriv_findiff2(a,dir,dx)
168 
169  s = numpy.shape(a) #size of the input array
170  aprime = numpy.array(a) #output has the same size than input
171 
172  #
173  # derivative for user dir=1, in python this is third coordinate
174  #
175  if dir == 1:
176  if (s[2] < 2):
177  return numpy.zeros_like(a)
178  if (s[2] < 4):
179  return deriv_findiff2(a,dir,dx)
180  if (s[2] < 6):
181  return deriv_findiff4(a,dir,dx)
182 
183  #forward differences near the first boundary
184  for i in range(3):
185  aprime[:,:,i] = (-11*a[:,:,i]+18*a[:,:,i+1]-9*a[:,:,i+2]+2*a[:,:,i+3]) / (6*dx)
186 
187  #centered differences
188  for i in range(3,s[2]-3):
189  aprime[:,:,i] = (-a[:,:,i-3]+9*a[:,:,i-2]-45*a[:,:,i-1]+45*a[:,:,i+1] -9*a[:,:,i+2]+a[:,:,i+3])/(60*dx)
190 
191  #backward differences near the second boundary
192  for i in range(s[2]-3,s[2]):
193  aprime[:,:,i] = (-2*a[:,:,i-3]+9*a[:,:,i-2]-18*a[:,:,i-1]+11*a[:,:,i]) /(6*dx)
194 
195  #
196  # derivative for dir=2
197  #
198  if dir == 2:
199  if (s[1] < 2):
200  return numpy.zeros_like(a)
201  if (s[1] < 4):
202  return deriv_findiff2(a,dir,dx)
203  if (s[1] < 6):
204  return deriv_findiff4(a,dir,dx)
205 
206  for i in range(3):
207  aprime[:,i,:] = (-11*a[:,i,:]+18*a[:,i+1,:]-9*a[:,i+2,:]+2*a[:,i+3,:]) /(6*dx) #forward differences near the first boundary
208 
209  for i in range(3,s[1]-3):
210  aprime[:,i,:] = (-a[:,i-3,:]+9*a[:,i-2,:]-45*a[:,i-1,:]+45*a[:,i+1,:] -9*a[:,i+2,:]+a[:,i+3,:])/(60*dx) #centered differences
211 
212  for i in range(s[1]-3,s[1]):
213  aprime[:,i,:] = (-2*a[:,i-3,:]+9*a[:,i-2,:]-18*a[:,i-1,:]+11*a[:,i,:]) /(6*dx) #backward differences near the second boundary
214 
215  #
216  # derivative for user dir=3
217  #
218  if dir == 3:
219  if (s[0] < 2):
220  return numpy.zeros_like(a)
221  if (s[0] < 4):
222  return deriv_findiff2(a,dir,dx)
223  if (s[0] < 6):
224  return deriv_findiff4(a,dir,dx)
225 
226  for i in range(3):
227  aprime[i,:,:] = (-11*a[i,:,:]+18*a[i+1,:,:]-9*a[i+2,:,:]+2*a[i+3,:,:]) /(6*dx) #forward differences near the first boundary
228 
229  for i in range(3,s[0]-3):
230  aprime[i,:,:] = (-a[i-3,:,:]+9*a[i-2,:,:]-45*a[i-1,:,:]+45*a[i+1,:,:] -9*a[i+2,:,:]+a[i+3,:,:])/(60*dx) #centered differences
231 
232  for i in range(s[0]-3,s[0]):
233  aprime[i,:,:] = (-2*a[i-3,:,:]+9*a[i-2,:,:]-18*a[i-1,:,:]+11*a[i,:,:]) /(6*dx) #backward differences near the second boundary
234 
235  return aprime
236 
237 def deriv_var_findiff2(a,var,dir):
238  ''' Function that calculates first-order derivatives
239  using 2nd order finite differences in regular Cartesian grids.'''
240 
241  import numpy
242  s = numpy.shape(a) #size of the input array
243  aprime = numpy.array(a) #output has the same size than input
244 
245  #
246  # derivative for user dir=1, in python this is third coordinate
247  #
248  if dir == 1:
249 
250  #forward differences near the first boundary
251  for i in range(1):
252  aprime[:,:,i] = (-a[:,:,i]+a[:,:,i+1]) / (-var[:,:,i]+var[:,:,i+1])
253 
254  #centered differences
255  for i in range(1,s[2]-1):
256  aprime[:,:,i] = (-a[:,:,i-1]+a[:,:,i+1])/(-var[:,:,i-1]+var[:,:,i+1])
257 
258  #backward differences near the second boundary
259  for i in range(s[2]-1,s[2]):
260  aprime[:,:,i] = (a[:,:,i-1]-a[:,:,i]) / (var[:,:,i-1]-var[:,:,i])
261 
262  #
263  # derivative for dir=2
264  #
265  if dir == 2:
266  #forward differences near the first boundary
267  for i in range(1):
268  aprime[:,i,:] = (-a[:,i,:]+a[:,i+1,:]) / (-var[:,i,:]+var[:,i+1,:])
269 
270  #centered differences
271  for i in range(1,s[1]-1):
272  aprime[:,i,:] = (-a[:,i-1,:]+a[:,i+1,:])/(-var[:,i-1,:]+var[:,i+1,:])
273 
274  #backward differences near the second boundary
275  for i in range(s[1]-1,s[1]):
276  aprime[:,i,:] = (a[:,i-1,:]-a[:,i,:]) / (var[:,i-1,:]-var[:,i,:])
277 
278  #
279  # derivative for user dir=3
280  #
281  if dir == 3:
282  #forward differences near the first boundary
283  for i in range(1):
284  aprime[i,:,:] = (-a[i,:,:]+a[i+1,:,:]) / (-var[i,:,:]+var[i+1,:,:])
285 
286  #centered differences
287  for i in range(1,s[0]-1):
288  aprime[i,:,:] = (-a[i-1,:,:]+a[i+1,:,:])/ (-var[i-1,:,:]+var[i+1,:,:])
289 
290  #backward differences near the second boundary
291  for i in range(s[0]-1,s[0]):
292  aprime[i,:,:] = (a[i-1,:,:]-a[i,:,:]) / (var[i-1,:,:]-var[i,:,:])
293 
294  return aprime
295 
296 
297 
298 def deriv_var_findiff4(a,var,dir):
299  ''' Function that calculates first-order derivatives of a variable
300  with respect to another variable 'var'
301  using 4th order finite differences in regular Cartesian grids.'''
302 
303  import numpy
304  s = numpy.shape(a) #size of the input array
305  aprime = numpy.array(a) #output has the same size than input
306 
307  #
308  # derivative for user dir=1, in python this is third coordinate
309  #
310  if dir == 1:
311 
312  #forward differences near the first boundary
313  for i in range(2):
314  aprime[:,:,i] = (-3*a[:,:,i]+4*a[:,:,i+1]-a[:,:,i+2]) / (-3*var[:,:,i]+4*var[:,:,i+1]-var[:,:,i+2])
315 
316  #centered differences
317  for i in range(2,s[2]-2):
318  aprime[:,:,i] = (a[:,:,i-2]-8*a[:,:,i-1]+8*a[:,:,i+1]-a[:,:,i+2])/(var[:,:,i-2]-8*var[:,:,i-1]+8*var[:,:,i+1]-var[:,:,i+2])
319 
320  #backward differences near the second boundary
321  for i in range(s[2]-2,s[2]):
322  aprime[:,:,i] = (a[:,:,i-2]-4*a[:,:,i-1]+3*a[:,:,i]) / (var[:,:,i-2]-4*var[:,:,i-1]+3*var[:,:,i])
323 
324  #
325  # derivative for dir=2
326  #
327  if dir == 2:
328  #forward differences near the first boundary
329  for i in range(2):
330  aprime[:,i,:] = (-3*a[:,i,:]+4*a[:,i+1,:]-a[:,i+2,:]) / (-3*var[:,i,:]+4*var[:,i+1,:]-var[:,i+2,:])
331 
332  #centered differences
333  for i in range(2,s[1]-2):
334  aprime[:,i,:] = (a[:,i-2,:]-8*a[:,i-1,:]+8*a[:,i+1,:]-a[:,i+2,:])/ (var[:,i-2,:]-8*var[:,i-1,:]+8*var[:,i+1,:]-var[:,i+2,:])
335 
336  #backward differences near the second boundary
337  for i in range(s[1]-2,s[1]):
338  aprime[:,i,:] = (a[:,i-2,:]-4*a[:,i-1,:]+3*a[:,i,:]) / (var[:,i-2,:]-4*var[:,i-1,:]+3*var[:,i,:])
339 
340  #
341  # derivative for user dir=3
342  #
343  if dir == 3:
344  #forward differences near the first boundary
345  for i in range(2):
346  aprime[i,:,:] = (-3*a[i,:,:]+4*a[i+1,:,:]-a[i+2,:,:]) / (-3*var[i,:,:]+4*var[i+1,:,:]-var[i+2,:,:])
347 
348  #centered differences
349  for i in range(2,s[0]-2):
350  aprime[i,:,:] = (a[i-2,:,:]-8*a[i-1,:,:]+8*a[i+1,:,:]-a[i+2,:,:])/ (var[i-2,:,:]-8*var[i-1,:,:]+8*var[i+1,:,:]-var[i+2,:,:])
351 
352  #backward differences near the second boundary
353  for i in range(s[0]-2,s[0]):
354  aprime[i,:,:] = (a[i-2,:,:]-4*a[i-1,:,:]+3*a[i,:,:]) / (var[i-2,:,:]-4*var[i-1,:,:]+3*var[i,:,:])
355 
356 
357  return aprime
358 
359 def deriv_var_findiff(a,var,dir,order=6):
360  ''' Function that calculates first-order derivatives of a variable
361  with respect to another variable 'var'
362  using sixth order finite differences in regular Cartesian grids.
363  The variable var should be increasing or decreasing in the
364  specified coordinate.
365  Calling sequence: aprime = deriv_var_findiff(a,var,dir,order)
366  a and var are 3-dimensional float32 numpy arrays
367  dir is 1, 2, or 3 (for X, Y, or Z directions in user coords.
368  order is the accuracy order, one of (2,4 or 6)
369  Note that these correspond to Z,Y,X directions in python coords.
370 
371  Returned array 'aprime' is the derivative of a wrt var. '''
372 
373  import numpy
374 
375  if order == 4:
376  return deriv_var_findiff4(a,var,dir)
377 
378  if order == 2:
379  return deriv_var_findiff2(a,var,dir)
380 
381  s = numpy.shape(a) #size of the input array
382  aprime = numpy.array(a) #output has the same size than input
383 
384  #
385  # derivative for dir=1
386  #
387  if dir == 1:
388  if (s[2] < 2):
389  return numpy.zeros_like(a)
390  if (s[2] < 4):
391  return deriv_var_findiff2(a,var,dir)
392  if (s[2] < 6):
393  return deriv_var_findiff4(a,var,dir)
394 
395  #forward differences near the first boundary
396  for i in range(3):
397  aprime[:,:,i] = (-11*a[:,:,i]+18*a[:,:,i+1]-9*a[:,:,i+2]+2*a[:,:,i+3]) / (-11*var[:,:,i]+18*var[:,:,i+1]-9*var[:,:,i+2]+2*var[:,:,i+3])
398 
399  #centered differences
400  for i in range(3,s[2]-3):
401  aprime[:,:,i] = (-a[:,:,i-3]+9*a[:,:,i-2]-45*a[:,:,i-1]+45*a[:,:,i+1] -9*a[:,:,i+2]+a[:,:,i+3])/(-var[:,:,i-3]+9*var[:,:,i-2]-45*var[:,:,i-1]+45*var[:,:,i+1] -9*var[:,:,i+2]+var[:,:,i+3])
402 
403  #backward differences near the second boundary
404  for i in range(s[2]-3,s[2]):
405  aprime[:,:,i] = (-2*a[:,:,i-3]+9*a[:,:,i-2]-18*a[:,:,i-1]+11*a[:,:,i]) /(-2*var[:,:,i-3]+9*var[:,:,i-2]-18*var[:,:,i-1]+11*var[:,:,i])
406 
407  #
408  # derivative for dir=2
409  #
410  if dir == 2:
411  if (s[1] < 2):
412  return numpy.zeros_like(a)
413  if (s[1] < 4):
414  return deriv_var_findiff2(a,var,dir)
415  if (s[1] < 6):
416  return deriv_var_findiff4(a,var,dir)
417 
418  for i in range(3):
419  aprime[:,i,:] = (-11*a[:,i,:]+18*a[:,i+1,:]-9*a[:,i+2,:]+2*a[:,i+3,:]) /(-11*var[:,i,:]+18*var[:,i+1,:]-9*var[:,i+2,:]+2*var[:,i+3,:]) #forward differences near the first boundary
420 
421  for i in range(3,s[1]-3):
422  aprime[:,i,:] = (-a[:,i-3,:]+9*a[:,i-2,:]-45*a[:,i-1,:]+45*a[:,i+1,:] -9*a[:,i+2,:]+a[:,i+3,:])/(-var[:,i-3,:]+9*var[:,i-2,:]-45*var[:,i-1,:]+45*var[:,i+1,:] -9*var[:,i+2,:]+var[:,i+3,:]) #centered differences
423 
424  for i in range(s[1]-3,s[1]):
425  aprime[:,i,:] = (-2*a[:,i-3,:]+9*a[:,i-2,:]-18*a[:,i-1,:]+11*a[:,i,:]) /(-2*var[:,i-3,:]+9*var[:,i-2,:]-18*var[:,i-1,:]+11*var[:,i,:]) #backward differences near the second boundary
426 
427  #
428  # derivative for dir=3
429  #
430  if dir == 3:
431  if (s[0] < 2):
432  return numpy.zeros_like(a)
433  if (s[0] < 4):
434  return deriv_var_findiff2(a,var,dir)
435  if (s[0] < 6):
436  return deriv_var_findiff4(a,var,dir)
437 
438  for i in range(3):
439  aprime[i,:,:] = (-11*a[i,:,:]+18*a[i+1,:,:]-9*a[i+2,:,:]+2*a[i+3,:,:]) /(-11*var[i,:,:]+18*var[i+1,:,:]-9*var[i+2,:,:]+2*var[i+3,:,:]) #forward differences near the first boundary
440 
441  for i in range(3,s[0]-3):
442  aprime[i,:,:] = (-a[i-3,:,:]+9*a[i-2,:,:]-45*a[i-1,:,:]+45*a[i+1,:,:] -9*a[i+2,:,:]+a[i+3,:,:])/(-var[i-3,:,:]+9*var[i-2,:,:]-45*var[i-1,:,:]+45*var[i+1,:,:] -9*var[i+2,:,:]+var[i+3,:,:]) #centered differences
443 
444  for i in range(s[0]-3,s[0]):
445  aprime[i,:,:] = (-2*a[i-3,:,:]+9*a[i-2,:,:]-18*a[i-1,:,:]+11*a[i,:,:]) /(-2*var[i-3,:,:]+9*var[i-2,:,:]-18*var[i-1,:,:]+11*var[i,:,:]) #backward differences near the second boundary
446 
447  return aprime
448 
449 def curl_findiff(A,B,C,order=6):
450  ''' Operator that calculates the curl of a vector field
451  using 6th order finite differences, on a regular Cartesian grid.
452  Calling sequence: curlfield = curl_findiff(A,B,C,order)
453  Where:
454  A,B,C are three 3-dimensional float32 arrays that define a
455  vector field.
456  order is the accuracy order (6,4, or 2)
457  curlfield is a 3-tuple of 3-dimensional float32 arrays that is
458  returned by this operator.'''
459  import vapor
460  ext = (vapor.BOUNDS[3]-vapor.BOUNDS[0], vapor.BOUNDS[4]-vapor.BOUNDS[1],vapor.BOUNDS[5]-vapor.BOUNDS[2])
461  usrmax = vapor.MapVoxToUser([vapor.BOUNDS[3],vapor.BOUNDS[4],vapor.BOUNDS[5]],vapor.REFINEMENT,vapor.LOD)
462  usrmin = vapor.MapVoxToUser([vapor.BOUNDS[0],vapor.BOUNDS[1],vapor.BOUNDS[2]],vapor.REFINEMENT,vapor.LOD)
463  dx = (usrmax[2]-usrmin[2])/ext[2] # grid delta is dx in user coord system
464  dy = (usrmax[1]-usrmin[1])/ext[1]
465  dz = (usrmax[0]-usrmin[0])/ext[0]
466 
467  aux1 = deriv_findiff(C,2,dy,order) #x component of the curl
468  aux2 = deriv_findiff(B,3,dz,order)
469  outx = aux1-aux2
470 
471  aux1 = deriv_findiff(A,3,dz,order) #y component of the curl
472  aux2 = deriv_findiff(C,1,dx,order)
473  outy = aux1-aux2
474 
475  aux1 = deriv_findiff(B,1,dx,order) #z component of the curl
476  aux2 = deriv_findiff(A,2,dy,order)
477  outz = aux1-aux2
478 
479  return outx, outy, outz #return results in user coordinate order
480 
481 
482 # Calculate divergence
483 def div_findiff(A,B,C,order=6):
484  ''' Operator that calculates the divergence of a vector field
485  using 6th order finite differences.
486  Calling sequence: DIV = div_findiff(A,B,C)
487  Where:
488  A, B, and C are 3-dimensional float32 arrays defining a vector field.
489  order is the accuracy order, one of (6,4, or 2)
490  Resulting DIV is a 3-dimensional float3d array consisting of
491  the divergence of the triple (A,B,C).'''
492  import vapor
493  ext = (vapor.BOUNDS[3]-vapor.BOUNDS[0], vapor.BOUNDS[4]-vapor.BOUNDS[1],vapor.BOUNDS[5]-vapor.BOUNDS[2])
494  usrmax = vapor.MapVoxToUser([vapor.BOUNDS[3],vapor.BOUNDS[4],vapor.BOUNDS[5]],vapor.REFINEMENT,vapor.LOD)
495  usrmin = vapor.MapVoxToUser([vapor.BOUNDS[0],vapor.BOUNDS[1],vapor.BOUNDS[2]],vapor.REFINEMENT,vapor.LOD)
496  dx = (usrmax[2]-usrmin[2])/ext[2] # grid delta-x in user coord system
497  dy = (usrmax[1]-usrmin[1])/ext[1]
498  dz = (usrmax[0]-usrmin[0])/ext[0]
499 # in User coords, A,B,C are x,y,z components
500  return deriv_findiff(C,3,dz,order)+deriv_findiff(B,2,dy,order)+deriv_findiff(A,1,dx,order)
501 
502 
503 def grad_findiff(A,order=6):
504  ''' Operator that calculates the gradient of a scalar field
505  using 6th order finite differences.
506  Calling sequence: GRD = grad_findiff(A)
507  Where:
508  A is a float32 array defining a scalar field.
509  order is the accuracy order, one of (6,4,or 2)
510  Result GRD is a triple of 3 3-dimensional float3d arrays consisting of
511  the gradient of A.'''
512  import vapor
513  ext = (vapor.BOUNDS[3]-vapor.BOUNDS[0], vapor.BOUNDS[4]-vapor.BOUNDS[1],vapor.BOUNDS[5]-vapor.BOUNDS[2]) #note that BOUNDS are in python coord order
514  usrmax = vapor.MapVoxToUser([vapor.BOUNDS[3],vapor.BOUNDS[4],vapor.BOUNDS[5]],vapor.REFINEMENT,vapor.LOD)
515  usrmin = vapor.MapVoxToUser([vapor.BOUNDS[0],vapor.BOUNDS[1],vapor.BOUNDS[2]],vapor.REFINEMENT,vapor.LOD)
516  dx = (usrmax[2]-usrmin[2])/ext[2] #delta in python z coordinate, user x
517  dy = (usrmax[1]-usrmin[1])/ext[1]
518  dz = (usrmax[0]-usrmin[0])/ext[0]# user z
519 
520  aux1 = deriv_findiff(A,1,dx,order) #x component of the gradient (in python system)
521  aux2 = deriv_findiff(A,2,dy,order)
522  aux3 = deriv_findiff(A,3,dz,order)
523 
524  return aux1,aux2,aux3 # return in user coordinate (x,y,z) order
525 
526 # Method that vertically interpolates one 3D variable to a level determined by
527 # another variable. The second variable (PR) is typically pressure.
528 # The second variable must decrease
529 # as a function of z (elevation). The returned value is a 2D variable having
530 # values interpolated to the surface defined by PR = val
531 # Sweep array from bottom to top
532 def interp3d(A,PR,val):
533  import numpy
534  s = numpy.shape(PR) #size of the input arrays
535  ss = [s[1],s[2]] # shape of 2d arrays
536  interpVal = numpy.empty(ss,numpy.float32)
537  ratio = numpy.zeros(ss,numpy.float32)
538 
539  # the LEVEL value is determine the lowest level where P<=val
540  LEVEL = numpy.empty(ss,numpy.int32)
541  LEVEL[:,:] = -1 #value where PR<=val has not been found
542  for K in range(s[0]):
543  #LEVNEED is true if this is first time PR<val.
544  LEVNEED = numpy.logical_and(numpy.less(LEVEL,0), numpy.less(PR[K,:,:] , val))
545  LEVEL[LEVNEED]=K
546  ratio[LEVNEED] = (val-PR[K,LEVNEED])/(PR[K-1,LEVNEED]-PR[K,LEVNEED])
547  interpVal[LEVNEED] = ratio[LEVNEED]*A[K,LEVNEED]+(1-ratio[LEVNEED])*A[K-1,LEVNEED]
548  LEVNEED = numpy.greater(LEVEL,0)
549  # Set unspecified values to value of A at top of data:
550  LEVNEED = numpy.less(LEVEL,0)
551  interpVal[LEVNEED] = A[s[0]-1,LEVNEED]
552  return interpVal
553 
554 def vector_rotate(angleRad, latDeg, u, v):
555  '''Rotate and scale vectors u,v for integration on
556  lon-lat grid.
557  Calling sequence:
558  rotfield=vector_rotate(angleRad, latDeg, u,v)
559  Where:
560  angleRad: 2D var, rotation from East in radians
561  latDeg: 2D var, latitude in degrees
562  u,v: 3D vars, x,y components of a vector field
563  rotfield is a 2-tuple of 3-dimensional float32 arrays,
564  representing rotation of u,v, returned by this operator.
565  '''
566  import numpy
567  import math
568  umod = numpy.cos(angleRad)*u + numpy.sin(angleRad)*v
569  vmod = -numpy.sin(angleRad)*u + numpy.cos(angleRad)*v
570  umod = umod/numpy.cos(latDeg*math.pi/180.)
571  return umod,vmod
572 
573 
vapor.MapVoxToUser
def MapVoxToUser(voxcoord, refinementLevel, lod)
Definition: vapor.py:20
vapor_utils.vector_rotate
def vector_rotate(angleRad, latDeg, u, v)
Definition: vapor_utils.py:554
vapor_utils.deriv_findiff
def deriv_findiff
Definition: vapor_utils.py:151
vapor_utils.deriv_findiff4
def deriv_findiff4(a, dir, dx)
Definition: vapor_utils.py:91
vapor_utils.div_findiff
def div_findiff
Definition: vapor_utils.py:483
vapor_utils.mag2d
def mag2d(a1, a2)
Definition: vapor_utils.py:23
vapor_utils.deriv_var_findiff2
def deriv_var_findiff2(a, var, dir)
Definition: vapor_utils.py:237
vapor_utils.grad_findiff
def grad_findiff
Definition: vapor_utils.py:503
vapor_utils.deriv_findiff2
def deriv_findiff2(a, dir, dx)
Definition: vapor_utils.py:32
vapor_utils.interp3d
def interp3d(A, PR, val)
Definition: vapor_utils.py:532
vapor_utils.curl_findiff
def curl_findiff
Definition: vapor_utils.py:449
vapor_utils.deriv_var_findiff4
def deriv_var_findiff4(a, var, dir)
Definition: vapor_utils.py:298
vapor_utils.deriv_var_findiff
def deriv_var_findiff
Definition: vapor_utils.py:359
vapor_utils.mag3d
def mag3d(a1, a2, a3)
Definition: vapor_utils.py:14
Generated by   doxygen 1.8.10