/[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 4209 - (hide annotations)
Mon Feb 18 06:00:56 2013 UTC (6 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 3460 byte(s)
Fixed joint example file, added ER Mapper plot script, updated content file.

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     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