-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplot_precipitation_climatology.py
106 lines (78 loc) · 3.29 KB
/
plot_precipitation_climatology.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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
import argparse
import numpy
import matplotlib.pyplot as plt
import iris
import iris.plot as iplt
import iris.coord_categorisation
iris.FUTURE.netcdf_promote = True
import cmocean
def read_data(fname, month):
"""Read an input data file"""
cube = iris.load_cube(fname, 'precipitation_flux')
iris.coord_categorisation.add_month(cube, 'time')
cube = cube.extract(iris.Constraint(month=month))
return cube
def convert_pr_units(cube):
"""Convert kg m-2 s-1 to mm day-1"""
expected_units = 'kg m-2 s-1'
assert cube.units == expected_units, 'Assumption is that units are {}'.format(expected_units)
cube.data = cube.data * 86400
cube.units = 'mm/day'
return cube
def plot_data(cube, month, gridlines=False, levels=None, mask=None):
"""Plot the data."""
import pdb
if not levels:
levels = numpy.arange(0, 10)
fig = plt.figure(figsize=[12,5])
iplt.contourf(cube, cmap=cmocean.cm.haline_r,
levels=levels,
extend='max')
plt.gca().coastlines()
if gridlines:
plt.gca().gridlines()
cbar = plt.colorbar()
cbar.set_label(str(cube.units))
title = '%s precipitation climatology (%s)' %(cube.attributes['model_id'], month)
plt.title(title)
def return_mask_from_land_surface_file(sftlf_file=None,mask='ocean'):
cube_land = iris.load_cube(sftlf_file,'land_area_fraction')
ocean_mask = numpy.where(cube_land.data < 5, True, False)
assert mask in ['ocean','land'],'mask must be ocean or land'
if mask == 'ocean':
return ocean_mask
elif mask == 'land':
return ~ocean_mask
def main(inargs):
"""Plot the precipitation climatology."""
cube = read_data(inargs.infile, inargs.month)
cube = convert_pr_units(cube)
clim = cube.collapsed('time', iris.analysis.MEAN)
if inargs.mask is not None:
clim.data = numpy.ma.asarray(clim.data)
clim.data.mask = return_mask_from_land_surface_file(
sftlf_file = inargs.mask[0],
mask = inargs.mask[1],
)
plot_data(clim, inargs.month, gridlines=inargs.gridlines,
levels=inargs.cbar_levels, mask=inargs.mask)
plt.savefig(inargs.outfile)
if __name__ == '__main__':
description='Plot the precipitation climatology for a given month.'
parser = argparse.ArgumentParser(description=description)
parser.add_argument("infile", type=str, help="Input file name")
parser.add_argument("month", type=str,
choices=['Jan', 'Feb', 'Mar', 'Apr',
'May', 'Jun', 'Jul', 'Aug',
'Sep', 'Oct', 'Nov', 'Dec'],
help="Month to plot")
parser.add_argument("outfile", type=str, help="Output file name")
parser.add_argument("--gridlines", action="store_true", default=False,
help="Include gridlines on the plot")
parser.add_argument("--cbar_levels", type=float, nargs='*', default=None,
help='list of levels / tick marks to appear on the colourbar')
parser.add_argument("--mask", type=str, nargs=2,
metavar=('sftlf_file', 'realm'), default=None,
help='Apply a land or ocean mask (specify the realm to mask)')
args = parser.parse_args()
main(args)