/[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 5593 - (show annotations)
Fri Apr 24 01:36:26 2015 UTC (4 years ago) by jfenwick
File MIME type: text/x-python
File size: 4100 byte(s)
Fixing institution name to comply with policy
1 from __future__ import division
2 from __future__ import print_function
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2015 by The University of Queensland
6 # 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 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 """This example show how to display ER Mapper raster data with matplotlib"""
19
20 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 from matplotlib import pyplot as plt
25 import numpy as np
26 import sys
27
28 # input filename
29 if len(sys.argv)>1:
30 FILENAME=sys.argv[1]
31 else:
32 FILENAME='data/QLDWestGravity.ers'
33
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 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 raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
72 except KeyError:
73 print("Cell type not specified. Assuming IEEE4ByteReal.")
74 datatype=np.float32
75
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 latitude=np.linspace(originY, originY-spacingY*NY, NY, endpoint=True)
109 DATA=np.fromfile(FILENAME, dtype=datatype).reshape(NY, NX)
110 # 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 plt.figure(figsize=(6*(spacingX*NX/(spacingY*NY))+1, 6), dpi=100)
115 plt.pcolor(x, y, DATA)
116 locs,_=plt.xticks()
117 plt.xticks(locs, list(map(lambda x:"%g"%x, locs)))
118 plt.xlabel(labelX)
119 plt.ylabel(labelY)
120 plt.axis('tight')
121 plt.title(FILENAME)
122 plt.colorbar()
123
124 plt.show()
125 plt.savefig("ermapper_plot.png")
126

  ViewVC Help
Powered by ViewVC 1.1.26