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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 4079 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2015 by The 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 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15
16 """This example shows how to display ER Mapper raster data with matplotlib"""
17 from __future__ import division, print_function
18
19 import matplotlib
20 # The following line is here to allow automated testing. Remove or comment if
21 # you would like to display the final plot in a window instead.
22 matplotlib.use('agg')
23 from matplotlib import pyplot as plt
24 import numpy as np
25 import sys
26
27 # input filename
28 if len(sys.argv)>1:
29 FILENAME=sys.argv[1]
30 else:
31 FILENAME='data/QLDWestGravity.ers'
32
33
34 if FILENAME[-4:]=='.ers': FILENAME=FILENAME[:-4]
35 ersfn=FILENAME+'.ers'
36 metadata=open(ersfn,'r').readlines()
37 # parse metadata
38 start=-1
39 for i in range(len(metadata)):
40 if metadata[i].strip() == 'DatasetHeader Begin':
41 start=i+1
42 if start==-1:
43 raise RuntimeError('Invalid ERS file (DatasetHeader not found)')
44
45 md_dict={}
46 section=[]
47 for i in range(start, len(metadata)):
48 line=metadata[i].strip()
49 if line[-6:].strip() == 'Begin':
50 section.append(line[:-6].strip())
51 elif line[-4:].strip() == 'End':
52 if len(section)>0:
53 section.pop()
54 else:
55 vals=line.split('=')
56 if len(vals)==2:
57 key = vals[0].strip()
58 value = vals[1].strip()
59 fullkey='.'.join(section+[key])
60 md_dict[fullkey]=value
61
62 try:
63 if md_dict['RasterInfo.CellType'] == 'IEEE4ByteReal':
64 datatype=np.float32
65 elif md_dict['RasterInfo.CellType'] == 'IEEE8ByteReal':
66 datatype=np.float64
67 elif md_dict['RasterInfo.CellType'] == 'Signed32BitInteger':
68 datatype=np.int32
69 else:
70 raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
71 except KeyError:
72 print("Cell type not specified. Assuming IEEE4ByteReal.")
73 datatype=np.float32
74
75 try:
76 NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
77 NY = int(md_dict['RasterInfo.NrOfLines'])
78 except:
79 raise RuntimeError("Could not determine extents of data")
80
81 try:
82 spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
83 spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
84 except:
85 raise RuntimeError("Could not determine cell dimensions")
86
87 try:
88 if md_dict['CoordinateSpace.CoordinateType']=='EN':
89 originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
90 originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
91 labelX = "Easting"
92 labelY = "Northing"
93 elif md_dict['CoordinateSpace.CoordinateType']=='LL':
94 originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
95 originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
96 labelX = "Longitude"
97 labelY = "Latitude"
98 else:
99 raise RuntimeError("Unknown CoordinateType")
100 except:
101 self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
102 originX,originY = 0.0, 0.0
103
104 f=open(FILENAME,'r')
105
106 longitude=np.linspace(originX, originX+spacingX*NX, NX, endpoint=True)
107 latitude=np.linspace(originY, originY-spacingY*NY, NY, endpoint=True)
108 DATA=np.fromfile(FILENAME, dtype=datatype).reshape(NY, NX)
109 # flip data in y-direction since ER Mapper stores data bottom up
110 DATA=np.flipud(DATA)
111
112 x,y=np.meshgrid(longitude, latitude)
113 plt.figure(figsize=(6*(spacingX*NX/(spacingY*NY))+1, 6), dpi=100)
114 plt.pcolor(x, y, DATA)
115 locs,_=plt.xticks()
116 plt.xticks(locs, list(map(lambda x:"%g"%x, locs)))
117 plt.xlabel(labelX)
118 plt.ylabel(labelY)
119 plt.axis('tight')
120 plt.title(FILENAME)
121 plt.colorbar()
122
123 plt.show()
124 plt.savefig("ermapper_plot.png")
125

  ViewVC Help
Powered by ViewVC 1.1.26