1 | #! /usr/bin/env python |
---|
2 | |
---|
3 | from ppclass import pp |
---|
4 | import numpy as np |
---|
5 | from scipy import ndimage |
---|
6 | |
---|
7 | dx = 50. ; filefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.nc" |
---|
8 | lim = 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 |
---|
16 | psfc = pp() |
---|
17 | psfc.file = filefile |
---|
18 | psfc.var = "PSFC" |
---|
19 | psfc.xcoeff = dx/1000. |
---|
20 | psfc.ycoeff = dx/1000. |
---|
21 | psfc.xlabel = "x distance (km)" |
---|
22 | psfc.ylabel = "y distance (km)" |
---|
23 | psfc.title = "Surface pressure $P_s$" |
---|
24 | #psfc.vmin = 482.5 |
---|
25 | #psfc.vmax = 483.5 |
---|
26 | psfc.div = 20 |
---|
27 | #psfc.out = "png" ; psfc.includedate = False ; psfc.folder = "dd/" |
---|
28 | |
---|
29 | ## calculate std |
---|
30 | std = np.std(psfc.getf()) |
---|
31 | |
---|
32 | ## calculate mean pressure |
---|
33 | ## ... same results whether or not time average |
---|
34 | psfcmean = pp() |
---|
35 | psfcmean << psfc |
---|
36 | psfcmean.x="0,10000" |
---|
37 | psfcmean.y="0,10000" |
---|
38 | psfcmean.t="0,10000" |
---|
39 | psfcmean.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 | |
---|
46 | for 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() |
---|