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

Annotation of /branches/doubleplusgood/doc/examples/inversion/plot_ermapper.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26