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

Annotation of /trunk/doc/examples/inversion/plot_ermapper.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4228 - (hide annotations)
Wed Feb 20 03:34:07 2013 UTC (6 years, 7 months ago) by caltinay
File MIME type: text/x-python
File size: 3528 byte(s)
Tweak & fix for data plotters

1 caltinay 4209 ##############################################################################
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 ER Mapper raster data with matplotlib"""
16    
17     from matplotlib import pyplot as plt
18     import numpy as np
19 caltinay 4228 import sys
20 caltinay 4209 from scipy.io import netcdf_file
21    
22     # input filename
23 caltinay 4228 if len(sys.argv)>1:
24     FILENAME=sys.argv[1]
25     else:
26     FILENAME='data/QLDWest_grav.ers'
27 caltinay 4209
28    
29     if FILENAME[-4:]=='.ers': FILENAME=FILENAME[:-4]
30     ersfn=FILENAME+'.ers'
31     metadata=open(ersfn,'r').readlines()
32     # parse metadata
33     start=-1
34     for i in range(len(metadata)):
35     if metadata[i].strip() == 'DatasetHeader Begin':
36     start=i+1
37     if start==-1:
38     raise RuntimeError('Invalid ERS file (DatasetHeader not found)')
39    
40     md_dict={}
41     section=[]
42     for i in range(start, len(metadata)):
43     line=metadata[i].strip()
44     if line[-6:].strip() == 'Begin':
45     section.append(line[:-6].strip())
46     elif line[-4:].strip() == 'End':
47     if len(section)>0:
48     section.pop()
49     else:
50     vals=line.split('=')
51     if len(vals)==2:
52     key = vals[0].strip()
53     value = vals[1].strip()
54     fullkey='.'.join(section+[key])
55     md_dict[fullkey]=value
56    
57     try:
58     if md_dict['RasterInfo.CellType'] != 'IEEE4ByteReal':
59     raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
60     except KeyError:
61     print("Cell type not specified. Assuming IEEE4ByteReal.")
62    
63     try:
64     NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
65     NY = int(md_dict['RasterInfo.NrOfLines'])
66     except:
67     raise RuntimeError("Could not determine extents of data")
68    
69     try:
70     spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
71     spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
72     except:
73     raise RuntimeError("Could not determine cell dimensions")
74    
75     try:
76     if md_dict['CoordinateSpace.CoordinateType']=='EN':
77     originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
78     originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
79     labelX = "Easting"
80     labelY = "Northing"
81     elif md_dict['CoordinateSpace.CoordinateType']=='LL':
82     originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
83     originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
84     labelX = "Longitude"
85     labelY = "Latitude"
86     else:
87     raise RuntimeError("Unknown CoordinateType")
88     except:
89     self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
90     originX,originY = 0.0, 0.0
91    
92     f=open(FILENAME,'r')
93    
94     longitude=np.linspace(originX, originX+spacingX*NX, NX, endpoint=True)
95 caltinay 4228 latitude=np.linspace(originY, originY-spacingY*NY, NY, endpoint=True)
96 caltinay 4209 DATA=np.fromfile(FILENAME, dtype=np.float32).reshape(NY, NX)
97     # flip data in y-direction since ER Mapper stores data bottom up
98     DATA=np.flipud(DATA)
99    
100     x,y=np.meshgrid(longitude, latitude)
101 caltinay 4228 plt.figure(figsize=(6*(spacingX*NX/(spacingY*NY))+1, 6), dpi=100)
102 caltinay 4209 plt.pcolor(x, y, DATA)
103     locs,_=plt.xticks()
104     plt.xticks(locs, map(lambda x:"%g"%x, locs))
105     plt.xlabel(labelX)
106     plt.ylabel(labelY)
107     plt.axis('tight')
108     plt.title(FILENAME)
109     plt.colorbar()
110     plt.show()
111    

  ViewVC Help
Powered by ViewVC 1.1.26