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

Last change on this file since 7672 was 7317, checked in by rwilson, 15 years ago

Replaced 'print' statements with log.critical() calls.

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