Changeset 1813 for inundation/pyvolution

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)

Location:

inundation/pyvolution

Files:

• ermapper_grids.py (modified) (5 diffs)
• test_ermapper.py (modified) (2 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    

inundation/pyvolution/test_ermapper.py

@@ -15 +15 @@
         pass
 
+    def test_write_grid(self):
+        header_filename = 'test_write_ermapper_grid.ers'
+        data_filename = 'test_write_ermapper_grid'
+        
+        original_grid = Numeric.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]])
+
+        # Check that the function works when passing the filename without
+        # a '.ers' extension
+        ermapper_grids.write_ermapper_grid(data_filename, original_grid)
+        new_grid = ermapper_grids.read_ermapper_grid(data_filename)
+        assert Numeric.allclose(original_grid, new_grid)
+
+        # Check that the function works when passing the filename with
+        # a '.ers' extension
+        ermapper_grids.write_ermapper_grid(header_filename, original_grid)
+        new_grid = ermapper_grids.read_ermapper_grid(header_filename)
+        assert Numeric.allclose(original_grid, new_grid)
+
+        # Clean up created files
+        remove(data_filename)
+        remove(header_filename)
+        
     def test_basic_single_line_grid(self):
         # Setup test data
@@ -56 +78 @@
         # Write test data
         ermapper_grids.write_ermapper_data(original_grid, data_filename)
-        # Write test header
+        # Write test header using all default values
         header_filename = data_filename + '.ers'                              
         ermapper_grids.write_ermapper_header(header_filename)
+
+        # Check that the read in values match the default values
+        header = ermapper_grids.read_ermapper_header(header_filename)
+
+        assert header['datum'] == '"GDA94"'
+        assert header['projection'] == '"GEOGRAPHIC"'  
+        assert header['coordinatetype'] == 'LL'
+        assert header['rotation'] == '0:0:0.0'
+        assert header['celltype'] == 'IEEE4ByteReal'
+        assert header['nullcellvalue'] == '-99999'
+        assert header['xdimension'] == '100'
+        assert header['registrationcellx'] == '0'
+        assert header['ydimension'] == '100'
+        assert header['registrationcelly'] == '2'
+        assert header['nroflines'] == '3'
+        assert header['nrofcellsperline'] == '4'
+        assert header['longitude'] == '0:0:0'
+        assert header['latitude'] == '0:0:0'
+        assert header['nrofbands'] == '1'
+        assert header['value'] == '"Default_Band"'
 
         # Clean up created files
         remove(data_filename)
         remove(header_filename)
+        
+    def test_header_creation(self):
+        header = {}
+        # have some values that aren't defaults
+        header['nroflines'] = '2'
+        header['nrofcellsperline'] = '3'
+
+        header = ermapper_grids.create_default_header(header)
+
+
+        # default values
+        assert header['datum'] == '"GDA94"'
+        assert header['projection'] == '"GEOGRAPHIC"'  
+        assert header['coordinatetype'] == 'LL'
+        assert header['rotation'] == '0:0:0.0'
+        assert header['celltype'] == 'IEEE4ByteReal'
+        assert header['nullcellvalue'] == '-99999'
+        assert header['xdimension'] == '100'
+        assert header['registrationcellx'] == '0'
+        assert header['ydimension'] == '100'
+        assert header['longitude'] == '0:0:0'
+        assert header['latitude'] == '0:0:0'
+        assert header['nrofbands'] == '1'
+        assert header['value'] == '"Default_Band"'
+
+        # non-default values
+        assert header['nroflines'] == '2'
+        assert header['nrofcellsperline'] == '3'
+        assert header['registrationcelly'] == '1'
+
+
+    def test_write_non_default_header(self):
+        header_filename = 'test_write_ermapper_grid.ers'
+
+        # setup test data
+        header = {}
+        # have some values that aren't defaults
+        header['nroflines'] = '2'
+        header['nrofcellsperline'] = '3'
+        
+        # Write test header using non-default values                           
+        ermapper_grids.write_ermapper_header(header_filename, header)
+
+        # Check that the read in values match the default values
+        header = ermapper_grids.read_ermapper_header(header_filename)
+
+        # default values
+        assert header['datum'] == '"GDA94"'
+        assert header['projection'] == '"GEOGRAPHIC"'  
+        assert header['coordinatetype'] == 'LL'
+        assert header['rotation'] == '0:0:0.0'
+        assert header['celltype'] == 'IEEE4ByteReal'
+        assert header['nullcellvalue'] == '-99999'
+        assert header['xdimension'] == '100'
+        assert header['registrationcellx'] == '0'
+        assert header['ydimension'] == '100'
+        assert header['longitude'] == '0:0:0'
+        assert header['latitude'] == '0:0:0'
+        assert header['nrofbands'] == '1'
+        assert header['value'] == '"Default_Band"'
+
+        # non-default values
+        assert header['nroflines'] == '2'
+        assert header['nrofcellsperline'] == '3'
+        assert header['registrationcelly'] == '1'
+
+        # Clean up created files
+        remove(header_filename)        
+
 
 # def test_default_filenames