/[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 4299 - (show annotations)
Mon Mar 11 00:40:09 2013 UTC (6 years, 1 month ago) by caltinay
File MIME type: text/x-python
File size: 3713 byte(s)
Exit with a warning in example files that require unavailable scipy.
bollide->hamster
added warnings to scons for scipy and gdal

1 ##############################################################################
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 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 from matplotlib import pyplot as plt
22 import numpy as np
23 import sys
24
25 # input filename
26 if len(sys.argv)>1:
27 FILENAME=sys.argv[1]
28 else:
29 FILENAME='data/QLDWestGravity.ers'
30
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 latitude=np.linspace(originY, originY-spacingY*NY, NY, endpoint=True)
99 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 plt.figure(figsize=(6*(spacingX*NX/(spacingY*NY))+1, 6), dpi=100)
105 plt.pcolor(x, y, DATA)
106 locs,_=plt.xticks()
107 plt.xticks(locs, map(lambda x:"%g"%x, locs))
108 plt.xlabel(labelX)
109 plt.ylabel(labelY)
110 plt.axis('tight')
111 plt.title(FILENAME)
112 plt.colorbar()
113
114 plt.show()
115 plt.savefig("ermapper_plot.png")
116

  ViewVC Help
Powered by ViewVC 1.1.26