/[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 4233 - (show annotations)
Thu Feb 21 02:00:09 2013 UTC (6 years, 7 months ago) by caltinay
File MIME type: text/x-python
File size: 1912 byte(s)
updated content file (wishlist). Added byte swapping to plotter.

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 from matplotlib import pyplot as plt
18 import numpy as np
19 import sys
20 from scipy.io import netcdf_file
21
22 # input filename
23 if len(sys.argv)>1:
24 FILENAME=sys.argv[1]
25 else:
26 FILENAME='data/QLDWest_grav.nc'
27
28 f=netcdf_file(FILENAME, 'r')
29 NY=f.dimensions["latitude"]
30 NX=f.dimensions["longitude"]
31
32 latitude=f.variables["latitude"]
33 y_label=latitude.long_name
34 y_units=latitude.units
35 latitude=latitude[:]
36
37 longitude=f.variables["longitude"]
38 x_label=longitude.long_name
39 x_units=longitude.units
40 longitude=longitude[:]
41
42 DATA=f.variables["onshore_only_Bouguer_geodetic"]
43 data_label=DATA.long_name
44 UNITS=DATA.units
45 DATA=DATA[:]
46 if DATA.dtype.byteorder=='>':
47 DATA=DATA.byteswap().newbyteorder()
48
49 f.close()
50
51 # add one more point in each dimension (because we fill cells):
52 ll=2*longitude[-1]-longitude[-2]
53 longitude=np.resize(longitude, len(longitude)+1)
54 longitude[-1]=ll
55 ll=2*latitude[-1]-latitude[-2]
56 latitude=np.resize(latitude, len(latitude)+1)
57 latitude[-1]=ll
58
59 lx=abs(longitude[-1]-longitude[0])
60 ly=abs(latitude[-1]-latitude[0])
61 x,y=np.meshgrid(longitude, latitude)
62 plt.figure(figsize=(6*lx/ly+1, 6), dpi=100)
63 plt.pcolor(x, y, DATA)
64 locs,_=plt.xticks()
65 plt.xticks(locs, map(lambda x:"%g"%x, locs))
66 plt.xlabel(x_label)
67 plt.ylabel(y_label)
68 plt.axis('tight')
69 plt.title("%s (%s)"%(data_label,UNITS))
70 plt.colorbar()
71 plt.show()
72

  ViewVC Help
Powered by ViewVC 1.1.26