/[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 4853 - (hide annotations)
Wed Apr 9 05:41:57 2014 UTC (5 years ago) by jfenwick
File MIME type: text/x-python
File size: 4096 byte(s)
Added import division and changed a few / to //

1 jfenwick 4853 from __future__ import division
2 jfenwick 4848 from __future__ import print_function
3 caltinay 4209 ##############################################################################
4     #
5 jfenwick 4657 # Copyright (c) 2003-2014 by University of Queensland
6 caltinay 4209 # http://www.uq.edu.au
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
14     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 caltinay 4209 #
16     ##############################################################################
17    
18     """This example show how to display ER Mapper raster data with matplotlib"""
19    
20 caltinay 4299 import matplotlib
21     # The following line is here to allow automated testing. Remove or comment if
22     # you would like to display the final plot in a window instead.
23     matplotlib.use('agg')
24 caltinay 4209 from matplotlib import pyplot as plt
25     import numpy as np
26 caltinay 4228 import sys
27 caltinay 4209
28     # input filename
29 caltinay 4228 if len(sys.argv)>1:
30     FILENAME=sys.argv[1]
31     else:
32 caltinay 4267 FILENAME='data/QLDWestGravity.ers'
33 caltinay 4209
34    
35     if FILENAME[-4:]=='.ers': FILENAME=FILENAME[:-4]
36     ersfn=FILENAME+'.ers'
37     metadata=open(ersfn,'r').readlines()
38     # parse metadata
39     start=-1
40     for i in range(len(metadata)):
41     if metadata[i].strip() == 'DatasetHeader Begin':
42     start=i+1
43     if start==-1:
44     raise RuntimeError('Invalid ERS file (DatasetHeader not found)')
45    
46     md_dict={}
47     section=[]
48     for i in range(start, len(metadata)):
49     line=metadata[i].strip()
50     if line[-6:].strip() == 'Begin':
51     section.append(line[:-6].strip())
52     elif line[-4:].strip() == 'End':
53     if len(section)>0:
54     section.pop()
55     else:
56     vals=line.split('=')
57     if len(vals)==2:
58     key = vals[0].strip()
59     value = vals[1].strip()
60     fullkey='.'.join(section+[key])
61     md_dict[fullkey]=value
62    
63     try:
64 caltinay 4495 if md_dict['RasterInfo.CellType'] == 'IEEE4ByteReal':
65     datatype=np.float32
66     elif md_dict['RasterInfo.CellType'] == 'IEEE8ByteReal':
67     datatype=np.float64
68     elif md_dict['RasterInfo.CellType'] == 'Signed32BitInteger':
69     datatype=np.int32
70     else:
71 caltinay 4209 raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
72     except KeyError:
73     print("Cell type not specified. Assuming IEEE4ByteReal.")
74 caltinay 4495 datatype=np.float32
75 caltinay 4209
76     try:
77     NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
78     NY = int(md_dict['RasterInfo.NrOfLines'])
79     except:
80     raise RuntimeError("Could not determine extents of data")
81    
82     try:
83     spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
84     spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
85     except:
86     raise RuntimeError("Could not determine cell dimensions")
87    
88     try:
89     if md_dict['CoordinateSpace.CoordinateType']=='EN':
90     originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
91     originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
92     labelX = "Easting"
93     labelY = "Northing"
94     elif md_dict['CoordinateSpace.CoordinateType']=='LL':
95     originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
96     originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
97     labelX = "Longitude"
98     labelY = "Latitude"
99     else:
100     raise RuntimeError("Unknown CoordinateType")
101     except:
102     self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
103     originX,originY = 0.0, 0.0
104    
105     f=open(FILENAME,'r')
106    
107     longitude=np.linspace(originX, originX+spacingX*NX, NX, endpoint=True)
108 caltinay 4228 latitude=np.linspace(originY, originY-spacingY*NY, NY, endpoint=True)
109 caltinay 4495 DATA=np.fromfile(FILENAME, dtype=datatype).reshape(NY, NX)
110 caltinay 4209 # flip data in y-direction since ER Mapper stores data bottom up
111     DATA=np.flipud(DATA)
112    
113     x,y=np.meshgrid(longitude, latitude)
114 caltinay 4228 plt.figure(figsize=(6*(spacingX*NX/(spacingY*NY))+1, 6), dpi=100)
115 caltinay 4209 plt.pcolor(x, y, DATA)
116     locs,_=plt.xticks()
117 caltinay 4306 plt.xticks(locs, list(map(lambda x:"%g"%x, locs)))
118 caltinay 4209 plt.xlabel(labelX)
119     plt.ylabel(labelY)
120     plt.axis('tight')
121     plt.title(FILENAME)
122     plt.colorbar()
123    
124 caltinay 4299 plt.show()
125     plt.savefig("ermapper_plot.png")
126 jfenwick 4289

  ViewVC Help
Powered by ViewVC 1.1.26