source: anuga_core/source/anuga/abstract_2d_finite_volumes/ermapper_grids.py @ 4936

Last change on this file since 4936 was 2054, checked in by ole, 19 years ago

Fixed unit tests

File size: 9.9 KB
Line 
1#!/usr/bin/env python
2
3# from os import open, write, read
4import Numeric
5
6celltype_map = {'IEEE4ByteReal': Numeric.Float32, 'IEEE8ByteReal': Numeric.Float64}
7
8
9def write_ermapper_grid(ofile, data, header = {}):
10    """
11    write_ermapper_grid(ofile, data, header = {}):
12   
13    Function to write a 2D Numeric array to an ERMapper grid.  There are a series of conventions adopted within
14    this code, specifically:
15    1)  The registration coordinate for the data is the SW (or lower-left) corner of the data
16    2)  The registration coordinates refer to cell centres
17    3)  The data is a 2D Numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
18        where N is the last line and M is the last column
19    4)  There has been no testng of the use of a rotated grid.  Best to keep data in an NS orientation
20   
21    Input Parameters:
22    ofile:      string - filename for output (note the output will consist of two files
23                ofile and ofile.ers.  Either of these can be entered into this function
24    data:       Numeric.array - 2D array containing the data to be output to the grid
25    header:     dictionary - contains spatial information about the grid, in particular:
26                    header['datum'] datum for the data ('"GDA94"')
27                    header['projection'] - either '"GEOGRAPHIC"' or '"PROJECTED"'
28                    header['coordinatetype'] - either 'EN' (for eastings/northings) or
29                                                      'LL' (for lat/long)
30                    header['rotation'] - rotation of grid ('0:0:0.0')
31                    header['celltype'] - data type for writing data ('IEEE4ByteReal')
32                    header['nullcellvalue'] - value for null cells ('-99999')
33                    header['xdimension'] - cell size in x-dir in units dictated by 'coordinatetype' ('100')
34                    header['registrationcellx'] == '0'
35                    header['ydimension'] - cell size in y-dir in units dictated by 'coordinatetype' ('100')
36                    header['longitude'] - co-ordinate of registration cell ('0:0:0')
37                    header['latitude'] - co-ordinate of registration line ('0:0:0')
38                    header['nrofbands'] - number of bands ('1')
39                    header['value'] - name of grid ('"Default_Band"')
40
41                    Some entries are determined automatically from the data
42                    header['nroflines'] - number of lines in data
43                    header['nrofcellsperline'] - number of columns in data
44                    header['registrationcelly'] == last line of data
45
46    Written by Trevor Dhu, Geoscience Australia 2005
47    """
48    # extract filenames for header and data files from ofile
49    ers_index = ofile.find('.ers')
50    if ers_index > 0:
51        data_file = ofile[0:ers_index]
52        header_file = ofile
53    else:
54        data_file = ofile
55        header_file = ofile + '.ers'
56
57
58    # Check that the data is a 2 dimensional array
59    data_size = Numeric.shape(data)
60    assert len(data_size) == 2
61   
62    header['nroflines'] = str(data_size[0])
63    header['nrofcellsperline'] = str(data_size[1])
64
65
66    header = create_default_header(header)
67    write_ermapper_header(header_file, header)
68    write_ermapper_data(data, data_file, data_format = header['celltype'])
69
70   
71def read_ermapper_grid(ifile):
72    ers_index = ifile.find('.ers')
73    if ers_index > 0:
74        data_file = ifile[0:ers_index]
75        header_file = ifile
76    else:
77        data_file = ifile
78        header_file = ifile + '.ers'
79       
80    header = read_ermapper_header(header_file)
81
82    nroflines = int(header['nroflines'])
83    nrofcellsperlines = int(header['nrofcellsperline'])
84    data = read_ermapper_data(data_file)
85    data = Numeric.reshape(data,(nroflines,nrofcellsperlines))
86    return data
87   
88
89def write_ermapper_header(ofile, header = {}):
90   
91    header = create_default_header(header)
92    # Determine if the dataset is in lats/longs or eastings/northings and set header parameters
93    # accordingly
94    if header['coordinatetype'] == 'LL':
95        X_Class = 'Longitude'
96        Y_Class = 'Latitude'
97    elif header['coordinatetype'] == 'EN':
98        X_Class = 'Eastings'
99        Y_Class = 'Northings'
100
101    # open the header file for writing to
102    fid = open(ofile,'wt')
103
104    # Begin writing the header
105    fid.write('DatasetHeader Begin\n')
106    fid.write('\tVersion\t\t= "6.4"\n')
107    fid.write('\tDatasetType\t= ERStorage\n')
108    fid.write('\tDataType\t= Raster\n')
109    fid.write('\tByteOrder\t= LSBFirst\n')
110
111    # Write the coordinate space information
112    fid.write('\tCoordinateSpace Begin\n')
113    fid.write('\t\tDatum\t\t\t = ' + header['datum'] + '\n')
114    fid.write('\t\tProjection\t\t = ' + header['projection'] + '\n')
115    fid.write('\t\tCoordinateType\t = ' + header['coordinatetype'] + '\n')
116    fid.write('\t\tRotation\t\t = ' + header['rotation'] + '\n')
117    fid.write('\t\tUnits\t\t = ' + header['units'] + '\n')
118    fid.write('\tCoordinateSpace End\n')
119
120    # Write the raster information
121    fid.write('\tRasterInfo Begin\n')
122    fid.write('\t\tCellType\t\t\t = ' + header['celltype'] + '\n')
123    fid.write('\t\tNullCellValue\t\t = ' + header['nullcellvalue'] + '\n')
124    fid.write('\t\tRegistrationCellX\t\t = ' + header['registrationcellx'] +'\n')
125    fid.write('\t\tRegistrationCellY\t\t = ' + header['registrationcelly'] +'\n')
126    # Write the cellsize information
127    fid.write('\t\tCellInfo Begin\n')
128    fid.write('\t\t\tXDimension\t\t\t = ' + header['xdimension'] + '\n')
129    fid.write('\t\t\tYDimension\t\t\t = ' + header['ydimension'] + '\n')
130    fid.write('\t\tCellInfo End\n')
131    # Continue with wrting the raster information
132    fid.write('\t\tNrOfLines\t\t\t = ' + header['nroflines'] + '\n')
133    fid.write('\t\tNrOfCellsPerLine\t = ' + header['nrofcellsperline'] + '\n')
134    # Write the registration coordinate information
135    fid.write('\t\tRegistrationCoord Begin\n')
136    ###print X_Class
137    fid.write('\t\t\t' + X_Class + '\t\t\t = ' + header[X_Class.lower()] + '\n')
138    fid.write('\t\t\t' + Y_Class + '\t\t\t = ' + header[Y_Class.lower()] + '\n')
139    fid.write('\t\tRegistrationCoord End\n')
140    # Continue with wrting the raster information
141    fid.write('\t\tNrOfBands\t\t\t = ' + header['nrofbands'] + '\n')
142    fid.write('\t\tBandID Begin\n')
143    fid.write('\t\t\tValue\t\t\t\t = ' + header['value'] + '\n')
144    fid.write('\t\tBandID End\n')
145    fid.write('\tRasterInfo End\n')
146    fid.write('DatasetHeader End\n')
147   
148    fid.close
149
150def read_ermapper_header(ifile):
151    # function for reading an ERMapper header from file
152    header = {}
153
154    fid = open(ifile,'rt')
155    header_string = fid.readlines()
156    fid.close()
157
158    for line in header_string:
159        if line.find('=') > 0:
160            tmp_string = line.strip().split('=')
161            header[tmp_string[0].strip().lower()]= tmp_string[1].strip()
162
163    return header                     
164
165def write_ermapper_data(grid, ofile, data_format = Numeric.Float32):
166
167
168    try:
169        data_format = celltype_map[data_format]
170    except:
171        pass
172   
173
174    #if isinstance(data_format, basestring):
175    #    #celltype_map is defined at top of code
176    #    if celltype_map.has_key(data_format):
177    #        data_format = celltype_map[data_format]
178    #    else:
179    #        msg = 'Format %s is not yet defined by celltype_map' %data_format
180    #        raise msg
181       
182   
183    # Convert the array to data_format (default format is Float32)
184    grid_as_float = grid.astype(data_format)
185
186    # Convert array to a string for writing to output file
187    output_string = grid_as_float.tostring()
188
189    # open output file in a binary format and write the output string
190    fid = open(ofile,'wb')
191    fid.write(output_string)
192    fid.close()
193
194
195def read_ermapper_data(ifile, data_format = Numeric.Float32):
196    # open input file in a binary format and read the input string
197    fid = open(ifile,'rb')
198    input_string = fid.read()
199    fid.close()
200
201    # convert input string to required format (Note default format is Numeric.Float32)
202    grid_as_float = Numeric.fromstring(input_string,data_format)
203    return grid_as_float
204
205def create_default_header(header = {}):
206    # fill any blanks in a header dictionary with default values
207    # input parameters:
208    # header:   a dictionary containing fields that are not meant
209    #           to be filled with default values
210   
211 
212    if not header.has_key('datum'):
213        header['datum'] = '"GDA94"'
214    if not header.has_key('projection'):
215        header['projection'] = '"GEOGRAPHIC"'       
216    if not header.has_key('coordinatetype'):
217        header['coordinatetype'] = 'LL'
218    if not header.has_key('rotation'):
219        header['rotation'] = '0:0:0.0'
220    if not header.has_key('units'):
221        header['units'] = '"METERS"'
222    if not header.has_key('celltype'):
223        header['celltype'] = 'IEEE4ByteReal'
224    if not header.has_key('nullcellvalue'):
225        header['nullcellvalue'] = '-99999'
226    if not header.has_key('xdimension'):
227        header['xdimension'] = '100'
228    if not header.has_key('latitude'):
229        header['latitude'] = '0:0:0'
230    if not header.has_key('longitude'):
231        header['longitude'] = '0:0:0'       
232    if not header.has_key('ydimension'):
233        header['ydimension'] = '100'
234    if not header.has_key('nroflines'):
235        header['nroflines'] = '3'   
236    if not header.has_key('nrofcellsperline'):
237        header['nrofcellsperline'] = '4'
238    if not header.has_key('registrationcellx'):
239        header['registrationcellx'] = '0'
240    if not header.has_key('registrationcelly'):
241        header['registrationcelly'] = str(int(header['nroflines'])-1)
242    if not header.has_key('nrofbands'):
243        header['nrofbands'] = '1'
244    if not header.has_key('value'):
245        header['value'] = '"Default_Band"'
246
247
248    return header   
249                         
250   
251
252
253   
Note: See TracBrowser for help on using the repository browser.