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

Contents of /trunk/doc/examples/inversion/plot_netcdf.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 3013 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2015 by The 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 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15
16 """This example shows how to display netCDF input data with matplotlib"""
17
18 from __future__ import division, print_function
19 import matplotlib
20 # The following line is here to allow automated testing. Remove or comment if
21 # you would like to display the final plot in a window instead.
22 matplotlib.use('agg')
23 from matplotlib import pyplot as plt
24 import numpy as np
25 import sys
26 try:
27 from scipy.io import netcdf_file
28
29 # input filename
30 if len(sys.argv)>1:
31 FILENAME=sys.argv[1]
32 else:
33 FILENAME='data/QLDWestMagnetic.nc'
34
35 f=netcdf_file(FILENAME, 'r')
36 for latvar in "latitude", "lat", "y":
37 try:
38 NY=f.dimensions[latvar]
39 break
40 except KeyError:
41 pass
42 for lonvar in "longitude", "lon", "x":
43 try:
44 NX=f.dimensions[lonvar]
45 break
46 except KeyError:
47 pass
48
49 latitude=f.variables[latvar]
50 y_label=latitude.long_name
51 y_units=latitude.units
52 latitude=latitude[:]
53
54 longitude=f.variables[lonvar]
55 x_label=longitude.long_name
56 x_units=longitude.units
57 longitude=longitude[:]
58
59 DATA=f.variables["magmap_V5_2010"]
60 data_label=DATA.long_name
61 try:
62 UNITS=DATA.units
63 except AttributeError:
64 UNITS="nT"
65
66 # convert NaN (not-a-number) values to -1000 so the plotting works
67 DATA=np.where(np.isnan(DATA[:]), -1000, DATA[:])
68
69 # make sure data is in native format
70 if DATA.dtype.byteorder != '=':
71 DATA=np.array(DATA, dtype=np.float32)
72
73 # add one more point in each dimension (because we fill cells):
74 ll=2*longitude[-1]-longitude[-2]
75 longitude=np.resize(longitude, len(longitude)+1)
76 longitude[-1]=ll
77 ll=2*latitude[-1]-latitude[-2]
78 latitude=np.resize(latitude, len(latitude)+1)
79 latitude[-1]=ll
80
81 lx=abs(longitude[-1]-longitude[0])
82 ly=abs(latitude[-1]-latitude[0])
83 x,y=np.meshgrid(longitude, latitude)
84 plt.figure(figsize=(6*lx/ly+1, 6), dpi=100)
85 plt.pcolor(x, y, DATA)
86 locs,_=plt.xticks()
87 plt.xticks(locs, list(map(lambda x:"%g"%x, locs)))
88 plt.xlabel(x_label)
89 plt.ylabel(y_label)
90 plt.axis('tight')
91 plt.title("%s (%s)"%(data_label,UNITS))
92 plt.colorbar()
93
94 plt.show()
95 plt.savefig("netcdf_plot.png")
96 # Now we can close the file
97 f.close()
98
99 except ImportError:
100 print("The scipy module was not found but is required to read netCDF files. Exiting...")
101

  ViewVC Help
Powered by ViewVC 1.1.26