/[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 4267 - (show annotations)
Fri Mar 1 04:50:54 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/x-python
File size: 3530 byte(s)
Updated inversion example scripts, added ER Mapper data, checked results
and corrected some issues. Cookbook does not reflect the changes yet.

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

  ViewVC Help
Powered by ViewVC 1.1.26