/[escript]/branches/doubleplusgood/doc/examples/inversion/plot_netcdf.py
ViewVC logotype

Contents of /branches/doubleplusgood/doc/examples/inversion/plot_netcdf.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4345 - (show annotations)
Fri Mar 29 07:09:41 2013 UTC (5 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 2632 byte(s)
Spelling fixes
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2013 by University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development since 2012 by School of Earth Sciences
12 #
13 ##############################################################################
14
15 """This example show how to display netCDF input data with matplotlib"""
16
17 import matplotlib
18 # The following line is here to allow automated testing. Remove or comment if
19 # you would like to display the final plot in a window instead.
20 matplotlib.use('agg')
21 from matplotlib import pyplot as plt
22 import numpy as np
23 import sys
24 try:
25 from scipy.io import netcdf_file
26
27 # input filename
28 if len(sys.argv)>1:
29 FILENAME=sys.argv[1]
30 else:
31 FILENAME='data/QLDWestMagnetic.nc'
32
33 f=netcdf_file(FILENAME, 'r')
34 NY=f.dimensions["latitude"]
35 NX=f.dimensions["longitude"]
36
37 latitude=f.variables["latitude"]
38 y_label=latitude.long_name
39 y_units=latitude.units
40 latitude=latitude[:]
41
42 longitude=f.variables["longitude"]
43 x_label=longitude.long_name
44 x_units=longitude.units
45 longitude=longitude[:]
46
47 DATA=f.variables["magmap_V5_2010"]
48 data_label=DATA.long_name
49 try:
50 UNITS=DATA.units
51 except AttributeError:
52 UNITS="nT"
53
54 # convert NaN (not-a-number) values to -1000 so the plotting works
55 DATA=np.where(np.isnan(DATA[:]), -1000, DATA[:])
56
57 # make sure data is in native format
58 if DATA.dtype.byteorder != '=':
59 DATA=np.array(DATA, dtype=np.float32)
60
61 f.close()
62
63 # add one more point in each dimension (because we fill cells):
64 ll=2*longitude[-1]-longitude[-2]
65 longitude=np.resize(longitude, len(longitude)+1)
66 longitude[-1]=ll
67 ll=2*latitude[-1]-latitude[-2]
68 latitude=np.resize(latitude, len(latitude)+1)
69 latitude[-1]=ll
70
71 lx=abs(longitude[-1]-longitude[0])
72 ly=abs(latitude[-1]-latitude[0])
73 x,y=np.meshgrid(longitude, latitude)
74 plt.figure(figsize=(6*lx/ly+1, 6), dpi=100)
75 plt.pcolor(x, y, DATA)
76 locs,_=plt.xticks()
77 plt.xticks(locs, list(map(lambda x:"%g"%x, locs)))
78 plt.xlabel(x_label)
79 plt.ylabel(y_label)
80 plt.axis('tight')
81 plt.title("%s (%s)"%(data_label,UNITS))
82 plt.colorbar()
83
84 plt.show()
85 plt.savefig("netcdf_plot.png")
86
87 except ImportError:
88 print("The scipy module was not found but is required to read netCDF files. Exiting...")
89

  ViewVC Help
Powered by ViewVC 1.1.26