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

Contents of /branches/doubleplusgood/doc/examples/inversion/plot_ermapper.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4209 - (show annotations)
Mon Feb 18 06:00:56 2013 UTC (6 years, 2 months ago) by caltinay
Original Path: trunk/doc/examples/inversion/plot_ermapper.py
File MIME type: text/x-python
File size: 3460 byte(s)
Fixed joint example file, added ER Mapper plot script, updated content file.

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

  ViewVC Help
Powered by ViewVC 1.1.26