/[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 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 2 months ago) by jfenwick
File MIME type: text/x-python
File size: 4026 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

  ViewVC Help
Powered by ViewVC 1.1.26