Changeset 1813 for inundation/pyvolution/ermapper_grids.py

Timestamp:

Sep 9, 2005, 9:26:34 AM (19 years ago)

Author:

tdhu

Message:

This is the first working version of the ermapper_girds.py that will write and read both ERMapper grids and binary data. Currently the reading of ERMapper grids only returns an array (not the spatial information that is read from the header file). Code also needs more checks on parameters such as data format (i.e. consistency between binary data format and header file)

File:

• inundation/pyvolution/ermapper_grids.py (modified) (5 diffs)

inundation/pyvolution/ermapper_grids.py

@@ -4 +4 @@
 import Numeric
 
+def write_ermapper_grid(ofile, data, header = {}):
+    # function write_ermapper_grid(ofile, data, header = {}):
+    #
+    # Function to write a 2D Numeric array to an ERMapper grid
+    #
+    # Input Parameters:
+    # ofile
 
-def write_ermapper_header(ofile, datum = 'AGD66', projection = 'GEOGRAPHIC',
-                          coordinate_type = 'LL', rotation = '0:0:0.0',
-                          cell_type = 'IEEE4ByteReal', null_value = '-99999',
-                          x_dim = '100', x_ll = '0:0:0', y_dim = '100', y_ll = '0:0:0',
-                          num_lines = '3', num_cells = '4', num_bands = '1',
-                          band_id = 'band1'):
+    # extract filenames for header and data files from ofile
+    ers_index = ofile.find('.ers')
+    if ers_index > 0:
+        data_file = ofile[0:ers_index]
+        header_file = ofile
+    else:
+        data_file = ofile
+        header_file = ofile + '.ers'
 
+
+    # Check that the data is a 2 dimensional array
+    data_size = Numeric.shape(data)
+    assert len(data_size) == 2
+    
+    header['nroflines'] = str(data_size[0])
+    header['nrofcellsperline'] = str(data_size[1])
+
+    write_ermapper_data(data,data_file)
+    write_ermapper_header(header_file,header)
+    
+def read_ermapper_grid(ifile):
+
+    ers_index = ifile.find('.ers')
+    if ers_index > 0:
+        data_file = ifile[0:ers_index]
+        header_file = ifile
+    else:
+        data_file = ifile
+        header_file = ifile + '.ers'
+        
+    header = read_ermapper_header(header_file)
+
+    nroflines = int(header['nroflines'])
+    nrofcellsperlines = int(header['nrofcellsperline'])
+    data = read_ermapper_data(data_file)
+    data = Numeric.reshape(data,(nroflines,nrofcellsperlines))
+    return data
+    
+
+def write_ermapper_header(ofile, header = {}):
+    
+    header = create_default_header(header)
     # Determine if the dataset is in lats/longs or eastings/northings and set header parameters
     # accordingly
-    if coordinate_type is 'LL':
+    if header['coordinatetype'] == 'LL':
         X_Class = 'Longitude'
         Y_Class = 'Latitude'
-    elif coordinate_type is 'EN':
+    elif header['coordinatetype'] == 'EN':
         X_Class = 'Eastings'
         Y_Class = 'Northings'
@@ -26 +68 @@
     # Begin writing the header
     fid.write('DatasetHeader Begin\n')
-    fid.write('\tVersion\t\t\t= "6.4"\n')
-    fid.write('\tDatasetType\t\t= ERStorage\n')
-    fid.write('\tDataType\t\t= Raster\n')
-    fid.write('\tByteOrder\t\t= LSBFirst\n')
+    fid.write('\tVersion\t\t= "6.4"\n')
+    fid.write('\tDatasetType\t= ERStorage\n')
+    fid.write('\tDataType\t= Raster\n')
+    fid.write('\tByteOrder\t= LSBFirst\n')
 
     # Write the coordinate space information
     fid.write('\tCoordinateSpace Begin\n')
-    fid.write('\t\tDatum\t\t\t = "' + datum + '"\n')
-    fid.write('\t\tProjection\t\t = "' + projection + '"\n')
-    fid.write('\t\tCoordinateType\t = ' + coordinate_type + '\n')
-    fid.write('\t\tRotation\t\t = ' + rotation + '\n')
+    fid.write('\t\tDatum\t\t\t = ' + header['datum'] + '\n')
+    fid.write('\t\tProjection\t\t = ' + header['projection'] + '\n')
+    fid.write('\t\tCoordinateType\t = ' + header['coordinatetype'] + '\n')
+    fid.write('\t\tRotation\t\t = ' + header['rotation'] + '\n')
     fid.write('\tCoordinateSpace End\n')
 
     # Write the raster information
     fid.write('\tRasterInfo Begin\n')
-    fid.write('\t\tCellType\t\t\t = ' + cell_type + '\n')
-    fid.write('\t\tNullCellValue\t\t = ' + null_value + '\n')
-    fid.write('\t\tRegistrationCellX\t\t = 0\n')
-    fid.write('\t\tRegistrationCellY\t\t = ' + str(int(num_lines) - 1) +'\n')
+    fid.write('\t\tCellType\t\t\t = ' + header['celltype'] + '\n')
+    fid.write('\t\tNullCellValue\t\t = ' + header['nullcellvalue'] + '\n')
+    fid.write('\t\tRegistrationCellX\t\t = ' + header['registrationcellx'] +'\n')
+    fid.write('\t\tRegistrationCellY\t\t = ' + header['registrationcelly'] +'\n')
     # Write the cellsize information
     fid.write('\t\tCellInfo Begin\n')
-    fid.write('\t\t\tXDimension\t\t\t = ' + x_dim + '\n')
-    fid.write('\t\t\tYDimension\t\t\t = ' + y_dim + '\n')
+    fid.write('\t\t\tXDimension\t\t\t = ' + header['xdimension'] + '\n')
+    fid.write('\t\t\tYDimension\t\t\t = ' + header['ydimension'] + '\n')
     fid.write('\t\tCellInfo End\n')
     # Continue with wrting the raster information
-    fid.write('\t\tNrOfLines\t\t\t = ' + num_lines + '\n')
-    fid.write('\t\tNrOfCellsPerLine\t = ' + num_cells + '\n')
+    fid.write('\t\tNrOfLines\t\t\t = ' + header['nroflines'] + '\n')
+    fid.write('\t\tNrOfCellsPerLine\t = ' + header['nrofcellsperline'] + '\n')
     # Write the registration coordinate information
     fid.write('\t\tRegistrationCoord Begin\n')
-    fid.write('\t\t\t' + X_Class + '\t\t\t = ' + x_ll + '\n')
-    fid.write('\t\t\t' + Y_Class + '\t\t\t = ' + y_ll + '\n')
+    fid.write('\t\t\t' + X_Class + '\t\t\t = ' + header[X_Class.lower()] + '\n')
+    fid.write('\t\t\t' + Y_Class + '\t\t\t = ' + header[Y_Class.lower()] + '\n')
     fid.write('\t\tRegistrationCoord End\n')
     # Continue with wrting the raster information
-    fid.write('\t\tNrOfBands\t\t\t = ' + num_bands + '\n')
+    fid.write('\t\tNrOfBands\t\t\t = ' + header['nrofbands'] + '\n')
     fid.write('\t\tBandID Begin\n')
-    fid.write('\t\t\tValue\t\t\t\t = "' + band_id + '"\n')
+    fid.write('\t\t\tValue\t\t\t\t = ' + header['value'] + '\n')
     fid.write('\t\tBandID End\n')
     fid.write('\tRasterInfo End\n')
@@ -67 +109 @@
     
     fid.close
-                          
+
+def read_ermapper_header(ifile):
+    # function for reading an ERMapper header from file
+    header = {}
+
+    fid = open(ifile,'rt')
+    header_string = fid.readlines()
+    fid.close()
+
+    for line in header_string:
+        if line.find('=') > 0:
+            tmp_string = line.strip().split('=')
+            header[tmp_string[0].strip().lower()]= tmp_string[1].strip()
+
+    return header                      
 
 def write_ermapper_data(grid, ofile, data_format = Numeric.Float32):
@@ -92 +148 @@
     return grid_as_float
 
-def create_default_header(header):
+def create_default_header(header = {}):
     # fill any blanks in a header dictionary with default values
     # input parameters:
@@ -98 +154 @@
     #           to be filled with default values
     
-##    projection = 'GEOGRAPHIC',
-##                          coordinate_type = 'LL', rotation = '0:0:0.0',
-##                          cell_type = 'IEEE4ByteReal', null_value = '-99999',
-##                          x_dim = '100', x_ll = '0:0:0', y_dim = '100', y_ll = '0:0:0',
-##                          num_lines = '3', num_cells = '4', num_bands = '1',
-##                          band_id = 'band1')
+ 
     if not header.has_key('datum'):
-        header['datum'] = 'GDA94'
+        header['datum'] = '"GDA94"'
     if not header.has_key('projection'):
-        header['datum'] = 'GEOGRAPHIC'        
+        header['projection'] = '"GEOGRAPHIC"'        
+    if not header.has_key('coordinatetype'):
+        header['coordinatetype'] = 'LL'
+    if not header.has_key('rotation'):
+        header['rotation'] = '0:0:0.0'
+    if not header.has_key('celltype'):
+        header['celltype'] = 'IEEE4ByteReal'
+    if not header.has_key('nullcellvalue'):
+        header['nullcellvalue'] = '-99999'
+    if not header.has_key('xdimension'):
+        header['xdimension'] = '100'
+    if not header.has_key('latitude'):
+        header['latitude'] = '0:0:0'
+    if not header.has_key('longitude'):
+        header['longitude'] = '0:0:0'        
+    if not header.has_key('ydimension'):
+        header['ydimension'] = '100'
+    if not header.has_key('nroflines'):
+        header['nroflines'] = '3'    
+    if not header.has_key('nrofcellsperline'):
+        header['nrofcellsperline'] = '4'
+    if not header.has_key('registrationcellx'):
+        header['registrationcellx'] = '0'
+    if not header.has_key('registrationcelly'):
+        header['registrationcelly'] = str(int(header['nroflines'])-1)
+    if not header.has_key('nrofbands'):
+        header['nrofbands'] = '1'
+    if not header.has_key('value'):
+        header['value'] = '"Default_Band"'
+
 
     return header