source: trunk/UTIL/PYTHON/powerlaw/showdd.py @ 1242

Last change on this file since 1242 was 1053, checked in by aslmd, 11 years ago

UTIL PYTHON. tools to study dust devils in LES.

  • Property svn:executable set to *
File size: 2.3 KB
Line 
1#! /usr/bin/env python
2
3from ppclass import pp
4import numpy as np
5from scipy import ndimage
6
7dx = 50. ; filefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.nc"
8lim = 4.
9#lim = -3.0
10#lim = -3.2
11
12#dx = 10. ; filefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.nc"
13
14## getplot pressure field
15## ... generic settings
16psfc = pp()
17psfc.file = filefile
18psfc.var = "PSFC"
19psfc.xcoeff = dx/1000.
20psfc.ycoeff = dx/1000.
21psfc.xlabel = "x distance (km)"
22psfc.ylabel = "y distance (km)"
23psfc.title = "Surface pressure $P_s$"
24#psfc.vmin = 482.5
25#psfc.vmax = 483.5
26psfc.div = 20
27#psfc.out = "png" ; psfc.includedate = False ; psfc.folder = "dd/"
28
29## calculate std
30std = np.std(psfc.getf())
31
32## calculate mean pressure
33## ... same results whether or not time average
34psfcmean = pp()
35psfcmean << psfc
36psfcmean.x="0,10000"
37psfcmean.y="0,10000"
38psfcmean.t="0,10000"
39psfcmean.getplot()
40
41#for t in [300]:
42#for t in range(0,500,50):
43#for t in [390,395,400,405,410,415]:
44#for t in [178,180,182]:
45
46for t in range(265,275,1):
47
48   ## plot pressure field
49   psfc.t = t
50   psfc.filename = "dd"+str(int(t))
51   psfc.getplot(extraplot=1)
52
53   ## calculate anomaly normalized by std
54   anomaly = (psfc - psfcmean) / std
55     
56   ## plot anomaly
57   anomaly.vmax = lim
58   anomaly.vmin = -4.0
59   anomaly.div = 32
60   anomaly.colorb = "gist_earth_r"
61   anomaly.colorb = "gist_stern"
62   anomaly.plotin = psfc
63   anomaly.title = "Anomaly $P_s - <P_s>^{xy}$"
64   anomaly.units = "$\sigma$"
65   #anomaly.plot()
66
67   ## sobel transform
68   field = anomaly.request[0][0][0][0][0][0].field
69   sx = ndimage.sobel(field, axis=0, mode='constant') 
70   sy = ndimage.sobel(field, axis=1, mode='constant')
71   tmp = np.hypot(sx, sy)
72
73   ## more prominent circles if we mutiply by field
74   tmp = tmp*field
75   #lim = -12.
76   #tmp[tmp>lim]=0.
77   #tmp[tmp<lim]=1.
78
79   #tmp = np.exp(-tmp)   
80
81   anomaly.request[0][0][0][0][0][0].field = tmp
82   anomaly.title = "Sobel transform" 
83   anomaly.units = "dimless" ; anomaly.vmax = None ; anomaly.vmin = None ; anomaly.div = 10
84   anomaly.colorb = "binary_r"
85   #anomaly.vmax = np.max(tmp)
86   anomaly.plot()
87
88   ### laplace transform
89   #anomaly.request[0][0][0][0][0][0].field=ndimage.laplace(anomaly.request[0][0][0][0][0][0].field)
90   #anomaly.title = "Sobel + Laplace transform"
91   #anomaly.units = "dimless" ; anomaly.vmax = None ; anomaly.vmin = None ; anomaly.div = 10
92   #anomaly.plot()
Note: See TracBrowser for help on using the repository browser.