/[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 4279 - (show annotations)
Wed Mar 6 05:10:44 2013 UTC (6 years, 1 month ago) by caltinay
File MIME type: text/x-python
File size: 3497 byte(s)
Update to API doco and separate section for multiple data sets.

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 from matplotlib import pyplot as plt
18 import numpy as np
19 import sys
20
21 # input filename
22 if len(sys.argv)>1:
23 FILENAME=sys.argv[1]
24 else:
25 FILENAME='data/QLDWestGravity.ers'
26
27
28 if FILENAME[-4:]=='.ers': FILENAME=FILENAME[:-4]
29 ersfn=FILENAME+'.ers'
30 metadata=open(ersfn,'r').readlines()
31 # parse metadata
32 start=-1
33 for i in range(len(metadata)):
34 if metadata[i].strip() == 'DatasetHeader Begin':
35 start=i+1
36 if start==-1:
37 raise RuntimeError('Invalid ERS file (DatasetHeader not found)')
38
39 md_dict={}
40 section=[]
41 for i in range(start, len(metadata)):
42 line=metadata[i].strip()
43 if line[-6:].strip() == 'Begin':
44 section.append(line[:-6].strip())
45 elif line[-4:].strip() == 'End':
46 if len(section)>0:
47 section.pop()
48 else:
49 vals=line.split('=')
50 if len(vals)==2:
51 key = vals[0].strip()
52 value = vals[1].strip()
53 fullkey='.'.join(section+[key])
54 md_dict[fullkey]=value
55
56 try:
57 if md_dict['RasterInfo.CellType'] != 'IEEE4ByteReal':
58 raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
59 except KeyError:
60 print("Cell type not specified. Assuming IEEE4ByteReal.")
61
62 try:
63 NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
64 NY = int(md_dict['RasterInfo.NrOfLines'])
65 except:
66 raise RuntimeError("Could not determine extents of data")
67
68 try:
69 spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
70 spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
71 except:
72 raise RuntimeError("Could not determine cell dimensions")
73
74 try:
75 if md_dict['CoordinateSpace.CoordinateType']=='EN':
76 originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
77 originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
78 labelX = "Easting"
79 labelY = "Northing"
80 elif md_dict['CoordinateSpace.CoordinateType']=='LL':
81 originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
82 originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
83 labelX = "Longitude"
84 labelY = "Latitude"
85 else:
86 raise RuntimeError("Unknown CoordinateType")
87 except:
88 self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
89 originX,originY = 0.0, 0.0
90
91 f=open(FILENAME,'r')
92
93 longitude=np.linspace(originX, originX+spacingX*NX, NX, endpoint=True)
94 latitude=np.linspace(originY, originY-spacingY*NY, NY, endpoint=True)
95 DATA=np.fromfile(FILENAME, dtype=np.float32).reshape(NY, NX)
96 # flip data in y-direction since ER Mapper stores data bottom up
97 DATA=np.flipud(DATA)
98
99 x,y=np.meshgrid(longitude, latitude)
100 plt.figure(figsize=(6*(spacingX*NX/(spacingY*NY))+1, 6), dpi=100)
101 plt.pcolor(x, y, DATA)
102 locs,_=plt.xticks()
103 plt.xticks(locs, map(lambda x:"%g"%x, locs))
104 plt.xlabel(labelX)
105 plt.ylabel(labelY)
106 plt.axis('tight')
107 plt.title(FILENAME)
108 plt.colorbar()
109 plt.show()
110

  ViewVC Help
Powered by ViewVC 1.1.26