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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4495 - (hide annotations)
Fri Jul 5 02:19:47 2013 UTC (5 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 3958 byte(s)
Added support for more input data types in ER Mapper files.
Tests will follow soon.

1 caltinay 4209 ##############################################################################
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 caltinay 4299 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 caltinay 4209 from matplotlib import pyplot as plt
22     import numpy as np
23 caltinay 4228 import sys
24 caltinay 4209
25     # input filename
26 caltinay 4228 if len(sys.argv)>1:
27     FILENAME=sys.argv[1]
28     else:
29 caltinay 4267 FILENAME='data/QLDWestGravity.ers'
30 caltinay 4209
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 caltinay 4495 if md_dict['RasterInfo.CellType'] == 'IEEE4ByteReal':
62     datatype=np.float32
63     elif md_dict['RasterInfo.CellType'] == 'IEEE8ByteReal':
64     datatype=np.float64
65     elif md_dict['RasterInfo.CellType'] == 'Signed32BitInteger':
66     datatype=np.int32
67     else:
68 caltinay 4209 raise RuntimeError('Unsupported data type '+md_dict['RasterInfo.CellType'])
69     except KeyError:
70     print("Cell type not specified. Assuming IEEE4ByteReal.")
71 caltinay 4495 datatype=np.float32
72 caltinay 4209
73     try:
74     NX = int(md_dict['RasterInfo.NrOfCellsPerLine'])
75     NY = int(md_dict['RasterInfo.NrOfLines'])
76     except:
77     raise RuntimeError("Could not determine extents of data")
78    
79     try:
80     spacingX = float(md_dict['RasterInfo.CellInfo.Xdimension'])
81     spacingY = float(md_dict['RasterInfo.CellInfo.Ydimension'])
82     except:
83     raise RuntimeError("Could not determine cell dimensions")
84    
85     try:
86     if md_dict['CoordinateSpace.CoordinateType']=='EN':
87     originX = float(md_dict['RasterInfo.RegistrationCoord.Eastings'])
88     originY = float(md_dict['RasterInfo.RegistrationCoord.Northings'])
89     labelX = "Easting"
90     labelY = "Northing"
91     elif md_dict['CoordinateSpace.CoordinateType']=='LL':
92     originX = float(md_dict['RasterInfo.RegistrationCoord.Longitude'])
93     originY = float(md_dict['RasterInfo.RegistrationCoord.Latitude'])
94     labelX = "Longitude"
95     labelY = "Latitude"
96     else:
97     raise RuntimeError("Unknown CoordinateType")
98     except:
99     self.logger.warn("Could not determine coordinate origin. Setting to (0.0, 0.0)")
100     originX,originY = 0.0, 0.0
101    
102     f=open(FILENAME,'r')
103    
104     longitude=np.linspace(originX, originX+spacingX*NX, NX, endpoint=True)
105 caltinay 4228 latitude=np.linspace(originY, originY-spacingY*NY, NY, endpoint=True)
106 caltinay 4495 DATA=np.fromfile(FILENAME, dtype=datatype).reshape(NY, NX)
107 caltinay 4209 # flip data in y-direction since ER Mapper stores data bottom up
108     DATA=np.flipud(DATA)
109    
110     x,y=np.meshgrid(longitude, latitude)
111 caltinay 4228 plt.figure(figsize=(6*(spacingX*NX/(spacingY*NY))+1, 6), dpi=100)
112 caltinay 4209 plt.pcolor(x, y, DATA)
113     locs,_=plt.xticks()
114 caltinay 4306 plt.xticks(locs, list(map(lambda x:"%g"%x, locs)))
115 caltinay 4209 plt.xlabel(labelX)
116     plt.ylabel(labelY)
117     plt.axis('tight')
118     plt.title(FILENAME)
119     plt.colorbar()
120    
121 caltinay 4299 plt.show()
122     plt.savefig("ermapper_plot.png")
123 jfenwick 4289

  ViewVC Help
Powered by ViewVC 1.1.26