/
eof_analysis_sst.py
88 lines (73 loc) · 2.93 KB
/
eof_analysis_sst.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#################################################################################################
# This Python script plots EOF analysis of sea-surface temperature observation
#
# Ji-Woo Lee, LLNL, July 2017
#################################################################################################
import cdms2 as cdms
import cdutil
import cdtime
from eofs.cdms import Eof
import vcs
#===========================================================================================================
# DATA
#-----------------------------------------------------------------------------------------------------------
# Open file ---
data_path = '/clim_obs/obs/ocn/mo/tos/UKMETOFFICE-HadISST-v1-1/130122_HadISST_sst.nc' ## Put your file here
f = cdms.open(data_path)
# Set time period ---
start_year = 1980
end_year = 2000
start_time = cdtime.comptime(start_year)
end_time = cdtime.comptime(end_year)
# Load variable ---
d = f('sst',time=(start_time,end_time),longitude=(0,360),latitude=(-90,90)) # Provide proper variable name
# Reomove annual cycle ---
d_anom = cdutil.ANNUALCYCLE.departures(d)
# EOF (take only first variance mode...) ---
solver = Eof(d_anom, weights='area')
eof = solver.eofsAsCovariance(neofs=1)
pc = solver.pcs(npcs=1, pcscaling=1) # pcscaling=1: scaled to unit variance
# (divided by the square-root of their eigenvalue)
frac = solver.varianceFraction()
# Sign control if needed ---
eof = eof * -1
pc = pc * -1
#===========================================================================================================
# Plot
#-----------------------------------------------------------------------------------------------------------
# Create canvas ---
canvas = vcs.init()
canvas.open()
template = canvas.createtemplate()
# Turn off no-needed information -- prevent overlap
template.blank(['title','mean','min','max','dataname','crdate','crtime',
'units','zvalue','tvalue','xunits','yunits','xname','yname'])
# Color setup ---
canvas.setcolormap('bl_to_darkred')
iso = canvas.createisofill()
levels = [-1, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
levels.insert(0,-1.e20)
levels.append(1.e20)
iso.levels = levels
colors = vcs.getcolors(levels,colors=range(16,240))
iso.fillareacolors = colors
iso.missing = 0 # Set missing value color as same as background
# Map projection ---
p = vcs.createprojection()
p.type = 'robinson'
iso.projection = p
# Plot ---
canvas.plot(eof[0],iso,template)
# Title ---
plot_title = vcs.createtext()
plot_title.x = .5
plot_title.y = .97
plot_title.height = 23
plot_title.halign = 'center'
plot_title.valign = 'top'
plot_title.color='black'
percentage = str(round(float(frac[0]*100.),1)) + '%' # % with one floating number
plot_title.string = 'EOF first mode, HadISST('+str(start_year)+'-'+str(end_year)+'), '+percentage
canvas.plot(plot_title)
# Save output as image file ---
canvas.png('eof_analysis_sst.png')