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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4344 by jfenwick, Wed Feb 27 03:42:40 2013 UTC revision 4345 by jfenwick, Fri Mar 29 07:09:41 2013 UTC
# Line 14  Line 14 
14    
15  """This example show how to display netCDF input data with matplotlib"""  """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  from matplotlib import pyplot as plt
22  import numpy as np  import numpy as np
23  import sys  import sys
24  from scipy.io import netcdf_file  try:
25        from scipy.io import netcdf_file
26    
27  # input filename      # input filename
28  if len(sys.argv)>1:      if len(sys.argv)>1:
29      FILENAME=sys.argv[1]          FILENAME=sys.argv[1]
30  else:      else:
31      FILENAME='data/QLDWest_grav.nc'          FILENAME='data/QLDWestMagnetic.nc'
32    
33  f=netcdf_file(FILENAME, 'r')      f=netcdf_file(FILENAME, 'r')
34  NY=f.dimensions["latitude"]      NY=f.dimensions["latitude"]
35  NX=f.dimensions["longitude"]      NX=f.dimensions["longitude"]
36    
37  latitude=f.variables["latitude"]      latitude=f.variables["latitude"]
38  y_label=latitude.long_name      y_label=latitude.long_name
39  y_units=latitude.units      y_units=latitude.units
40  latitude=latitude[:]      latitude=latitude[:]
41    
42  longitude=f.variables["longitude"]      longitude=f.variables["longitude"]
43  x_label=longitude.long_name      x_label=longitude.long_name
44  x_units=longitude.units      x_units=longitude.units
45  longitude=longitude[:]      longitude=longitude[:]
46    
47  DATA=f.variables["onshore_only_Bouguer_geodetic"]      DATA=f.variables["magmap_V5_2010"]
48  data_label=DATA.long_name      data_label=DATA.long_name
49  UNITS=DATA.units      try:
50  DATA=DATA[:]          UNITS=DATA.units
51  if DATA.dtype.byteorder=='>':      except AttributeError:
52      DATA=DATA.byteswap().newbyteorder()          UNITS="nT"
53    
54  f.close()      # convert NaN (not-a-number) values to -1000 so the plotting works
55        DATA=np.where(np.isnan(DATA[:]), -1000, DATA[:])
56  # add one more point in each dimension (because we fill cells):  
57  ll=2*longitude[-1]-longitude[-2]      # make sure data is in native format
58  longitude=np.resize(longitude, len(longitude)+1)      if DATA.dtype.byteorder != '=':
59  longitude[-1]=ll          DATA=np.array(DATA, dtype=np.float32)
60  ll=2*latitude[-1]-latitude[-2]  
61  latitude=np.resize(latitude, len(latitude)+1)      f.close()
62  latitude[-1]=ll  
63        # add one more point in each dimension (because we fill cells):
64  lx=abs(longitude[-1]-longitude[0])      ll=2*longitude[-1]-longitude[-2]
65  ly=abs(latitude[-1]-latitude[0])      longitude=np.resize(longitude, len(longitude)+1)
66  x,y=np.meshgrid(longitude, latitude)      longitude[-1]=ll
67  plt.figure(figsize=(6*lx/ly+1, 6), dpi=100)      ll=2*latitude[-1]-latitude[-2]
68  plt.pcolor(x, y, DATA)      latitude=np.resize(latitude, len(latitude)+1)
69  locs,_=plt.xticks()      latitude[-1]=ll
70  plt.xticks(locs, map(lambda x:"%g"%x, locs))  
71  plt.xlabel(x_label)      lx=abs(longitude[-1]-longitude[0])
72  plt.ylabel(y_label)      ly=abs(latitude[-1]-latitude[0])
73  plt.axis('tight')      x,y=np.meshgrid(longitude, latitude)
74  plt.title("%s (%s)"%(data_label,UNITS))      plt.figure(figsize=(6*lx/ly+1, 6), dpi=100)
75  plt.colorbar()      plt.pcolor(x, y, DATA)
76  plt.show()      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    

Legend:
Removed from v.4344  
changed lines
  Added in v.4345

  ViewVC Help
Powered by ViewVC 1.1.26