Changeset 1813


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/ermapper_grids.py

    r1812 r1813  
    44import Numeric
    55
     6def write_ermapper_grid(ofile, data, header = {}):
     7    # function write_ermapper_grid(ofile, data, header = {}):
     8    #
     9    # Function to write a 2D Numeric array to an ERMapper grid
     10    #
     11    # Input Parameters:
     12    # ofile
    613
    7 def write_ermapper_header(ofile, datum = 'AGD66', projection = 'GEOGRAPHIC',
    8                           coordinate_type = 'LL', rotation = '0:0:0.0',
    9                           cell_type = 'IEEE4ByteReal', null_value = '-99999',
    10                           x_dim = '100', x_ll = '0:0:0', y_dim = '100', y_ll = '0:0:0',
    11                           num_lines = '3', num_cells = '4', num_bands = '1',
    12                           band_id = 'band1'):
     14    # extract filenames for header and data files from ofile
     15    ers_index = ofile.find('.ers')
     16    if ers_index > 0:
     17        data_file = ofile[0:ers_index]
     18        header_file = ofile
     19    else:
     20        data_file = ofile
     21        header_file = ofile + '.ers'
    1322
     23
     24    # Check that the data is a 2 dimensional array
     25    data_size = Numeric.shape(data)
     26    assert len(data_size) == 2
     27   
     28    header['nroflines'] = str(data_size[0])
     29    header['nrofcellsperline'] = str(data_size[1])
     30
     31    write_ermapper_data(data,data_file)
     32    write_ermapper_header(header_file,header)
     33   
     34def read_ermapper_grid(ifile):
     35
     36    ers_index = ifile.find('.ers')
     37    if ers_index > 0:
     38        data_file = ifile[0:ers_index]
     39        header_file = ifile
     40    else:
     41        data_file = ifile
     42        header_file = ifile + '.ers'
     43       
     44    header = read_ermapper_header(header_file)
     45
     46    nroflines = int(header['nroflines'])
     47    nrofcellsperlines = int(header['nrofcellsperline'])
     48    data = read_ermapper_data(data_file)
     49    data = Numeric.reshape(data,(nroflines,nrofcellsperlines))
     50    return data
     51   
     52
     53def write_ermapper_header(ofile, header = {}):
     54   
     55    header = create_default_header(header)
    1456    # Determine if the dataset is in lats/longs or eastings/northings and set header parameters
    1557    # accordingly
    16     if coordinate_type is 'LL':
     58    if header['coordinatetype'] == 'LL':
    1759        X_Class = 'Longitude'
    1860        Y_Class = 'Latitude'
    19     elif coordinate_type is 'EN':
     61    elif header['coordinatetype'] == 'EN':
    2062        X_Class = 'Eastings'
    2163        Y_Class = 'Northings'
     
    2668    # Begin writing the header
    2769    fid.write('DatasetHeader Begin\n')
    28     fid.write('\tVersion\t\t\t= "6.4"\n')
    29     fid.write('\tDatasetType\t\t= ERStorage\n')
    30     fid.write('\tDataType\t\t= Raster\n')
    31     fid.write('\tByteOrder\t\t= LSBFirst\n')
     70    fid.write('\tVersion\t\t= "6.4"\n')
     71    fid.write('\tDatasetType\t= ERStorage\n')
     72    fid.write('\tDataType\t= Raster\n')
     73    fid.write('\tByteOrder\t= LSBFirst\n')
    3274
    3375    # Write the coordinate space information
    3476    fid.write('\tCoordinateSpace Begin\n')
    35     fid.write('\t\tDatum\t\t\t = "' + datum + '"\n')
    36     fid.write('\t\tProjection\t\t = "' + projection + '"\n')
    37     fid.write('\t\tCoordinateType\t = ' + coordinate_type + '\n')
    38     fid.write('\t\tRotation\t\t = ' + rotation + '\n')
     77    fid.write('\t\tDatum\t\t\t = ' + header['datum'] + '\n')
     78    fid.write('\t\tProjection\t\t = ' + header['projection'] + '\n')
     79    fid.write('\t\tCoordinateType\t = ' + header['coordinatetype'] + '\n')
     80    fid.write('\t\tRotation\t\t = ' + header['rotation'] + '\n')
    3981    fid.write('\tCoordinateSpace End\n')
    4082
    4183    # Write the raster information
    4284    fid.write('\tRasterInfo Begin\n')
    43     fid.write('\t\tCellType\t\t\t = ' + cell_type + '\n')
    44     fid.write('\t\tNullCellValue\t\t = ' + null_value + '\n')
    45     fid.write('\t\tRegistrationCellX\t\t = 0\n')
    46     fid.write('\t\tRegistrationCellY\t\t = ' + str(int(num_lines) - 1) +'\n')
     85    fid.write('\t\tCellType\t\t\t = ' + header['celltype'] + '\n')
     86    fid.write('\t\tNullCellValue\t\t = ' + header['nullcellvalue'] + '\n')
     87    fid.write('\t\tRegistrationCellX\t\t = ' + header['registrationcellx'] +'\n')
     88    fid.write('\t\tRegistrationCellY\t\t = ' + header['registrationcelly'] +'\n')
    4789    # Write the cellsize information
    4890    fid.write('\t\tCellInfo Begin\n')
    49     fid.write('\t\t\tXDimension\t\t\t = ' + x_dim + '\n')
    50     fid.write('\t\t\tYDimension\t\t\t = ' + y_dim + '\n')
     91    fid.write('\t\t\tXDimension\t\t\t = ' + header['xdimension'] + '\n')
     92    fid.write('\t\t\tYDimension\t\t\t = ' + header['ydimension'] + '\n')
    5193    fid.write('\t\tCellInfo End\n')
    5294    # Continue with wrting the raster information
    53     fid.write('\t\tNrOfLines\t\t\t = ' + num_lines + '\n')
    54     fid.write('\t\tNrOfCellsPerLine\t = ' + num_cells + '\n')
     95    fid.write('\t\tNrOfLines\t\t\t = ' + header['nroflines'] + '\n')
     96    fid.write('\t\tNrOfCellsPerLine\t = ' + header['nrofcellsperline'] + '\n')
    5597    # Write the registration coordinate information
    5698    fid.write('\t\tRegistrationCoord Begin\n')
    57     fid.write('\t\t\t' + X_Class + '\t\t\t = ' + x_ll + '\n')
    58     fid.write('\t\t\t' + Y_Class + '\t\t\t = ' + y_ll + '\n')
     99    fid.write('\t\t\t' + X_Class + '\t\t\t = ' + header[X_Class.lower()] + '\n')
     100    fid.write('\t\t\t' + Y_Class + '\t\t\t = ' + header[Y_Class.lower()] + '\n')
    59101    fid.write('\t\tRegistrationCoord End\n')
    60102    # Continue with wrting the raster information
    61     fid.write('\t\tNrOfBands\t\t\t = ' + num_bands + '\n')
     103    fid.write('\t\tNrOfBands\t\t\t = ' + header['nrofbands'] + '\n')
    62104    fid.write('\t\tBandID Begin\n')
    63     fid.write('\t\t\tValue\t\t\t\t = "' + band_id + '"\n')
     105    fid.write('\t\t\tValue\t\t\t\t = ' + header['value'] + '\n')
    64106    fid.write('\t\tBandID End\n')
    65107    fid.write('\tRasterInfo End\n')
     
    67109   
    68110    fid.close
    69                          
     111
     112def read_ermapper_header(ifile):
     113    # function for reading an ERMapper header from file
     114    header = {}
     115
     116    fid = open(ifile,'rt')
     117    header_string = fid.readlines()
     118    fid.close()
     119
     120    for line in header_string:
     121        if line.find('=') > 0:
     122            tmp_string = line.strip().split('=')
     123            header[tmp_string[0].strip().lower()]= tmp_string[1].strip()
     124
     125    return header                     
    70126
    71127def write_ermapper_data(grid, ofile, data_format = Numeric.Float32):
     
    92148    return grid_as_float
    93149
    94 def create_default_header(header):
     150def create_default_header(header = {}):
    95151    # fill any blanks in a header dictionary with default values
    96152    # input parameters:
     
    98154    #           to be filled with default values
    99155   
    100 ##    projection = 'GEOGRAPHIC',
    101 ##                          coordinate_type = 'LL', rotation = '0:0:0.0',
    102 ##                          cell_type = 'IEEE4ByteReal', null_value = '-99999',
    103 ##                          x_dim = '100', x_ll = '0:0:0', y_dim = '100', y_ll = '0:0:0',
    104 ##                          num_lines = '3', num_cells = '4', num_bands = '1',
    105 ##                          band_id = 'band1')
     156 
    106157    if not header.has_key('datum'):
    107         header['datum'] = 'GDA94'
     158        header['datum'] = '"GDA94"'
    108159    if not header.has_key('projection'):
    109         header['datum'] = 'GEOGRAPHIC'       
     160        header['projection'] = '"GEOGRAPHIC"'       
     161    if not header.has_key('coordinatetype'):
     162        header['coordinatetype'] = 'LL'
     163    if not header.has_key('rotation'):
     164        header['rotation'] = '0:0:0.0'
     165    if not header.has_key('celltype'):
     166        header['celltype'] = 'IEEE4ByteReal'
     167    if not header.has_key('nullcellvalue'):
     168        header['nullcellvalue'] = '-99999'
     169    if not header.has_key('xdimension'):
     170        header['xdimension'] = '100'
     171    if not header.has_key('latitude'):
     172        header['latitude'] = '0:0:0'
     173    if not header.has_key('longitude'):
     174        header['longitude'] = '0:0:0'       
     175    if not header.has_key('ydimension'):
     176        header['ydimension'] = '100'
     177    if not header.has_key('nroflines'):
     178        header['nroflines'] = '3'   
     179    if not header.has_key('nrofcellsperline'):
     180        header['nrofcellsperline'] = '4'
     181    if not header.has_key('registrationcellx'):
     182        header['registrationcellx'] = '0'
     183    if not header.has_key('registrationcelly'):
     184        header['registrationcelly'] = str(int(header['nroflines'])-1)
     185    if not header.has_key('nrofbands'):
     186        header['nrofbands'] = '1'
     187    if not header.has_key('value'):
     188        header['value'] = '"Default_Band"'
     189
    110190
    111191    return header   
  • inundation/pyvolution/test_ermapper.py

    r1802 r1813  
    1515        pass
    1616
     17    def test_write_grid(self):
     18        header_filename = 'test_write_ermapper_grid.ers'
     19        data_filename = 'test_write_ermapper_grid'
     20       
     21        original_grid = Numeric.array([[0.0, 0.1, 1.0], [2.0, 3.0, 4.0]])
     22
     23        # Check that the function works when passing the filename without
     24        # a '.ers' extension
     25        ermapper_grids.write_ermapper_grid(data_filename, original_grid)
     26        new_grid = ermapper_grids.read_ermapper_grid(data_filename)
     27        assert Numeric.allclose(original_grid, new_grid)
     28
     29        # Check that the function works when passing the filename with
     30        # a '.ers' extension
     31        ermapper_grids.write_ermapper_grid(header_filename, original_grid)
     32        new_grid = ermapper_grids.read_ermapper_grid(header_filename)
     33        assert Numeric.allclose(original_grid, new_grid)
     34
     35        # Clean up created files
     36        remove(data_filename)
     37        remove(header_filename)
     38       
    1739    def test_basic_single_line_grid(self):
    1840        # Setup test data
     
    5678        # Write test data
    5779        ermapper_grids.write_ermapper_data(original_grid, data_filename)
    58         # Write test header
     80        # Write test header using all default values
    5981        header_filename = data_filename + '.ers'                             
    6082        ermapper_grids.write_ermapper_header(header_filename)
     83
     84        # Check that the read in values match the default values
     85        header = ermapper_grids.read_ermapper_header(header_filename)
     86
     87        assert header['datum'] == '"GDA94"'
     88        assert header['projection'] == '"GEOGRAPHIC"' 
     89        assert header['coordinatetype'] == 'LL'
     90        assert header['rotation'] == '0:0:0.0'
     91        assert header['celltype'] == 'IEEE4ByteReal'
     92        assert header['nullcellvalue'] == '-99999'
     93        assert header['xdimension'] == '100'
     94        assert header['registrationcellx'] == '0'
     95        assert header['ydimension'] == '100'
     96        assert header['registrationcelly'] == '2'
     97        assert header['nroflines'] == '3'
     98        assert header['nrofcellsperline'] == '4'
     99        assert header['longitude'] == '0:0:0'
     100        assert header['latitude'] == '0:0:0'
     101        assert header['nrofbands'] == '1'
     102        assert header['value'] == '"Default_Band"'
    61103
    62104        # Clean up created files
    63105        remove(data_filename)
    64106        remove(header_filename)
     107       
     108    def test_header_creation(self):
     109        header = {}
     110        # have some values that aren't defaults
     111        header['nroflines'] = '2'
     112        header['nrofcellsperline'] = '3'
     113
     114        header = ermapper_grids.create_default_header(header)
     115
     116
     117        # default values
     118        assert header['datum'] == '"GDA94"'
     119        assert header['projection'] == '"GEOGRAPHIC"' 
     120        assert header['coordinatetype'] == 'LL'
     121        assert header['rotation'] == '0:0:0.0'
     122        assert header['celltype'] == 'IEEE4ByteReal'
     123        assert header['nullcellvalue'] == '-99999'
     124        assert header['xdimension'] == '100'
     125        assert header['registrationcellx'] == '0'
     126        assert header['ydimension'] == '100'
     127        assert header['longitude'] == '0:0:0'
     128        assert header['latitude'] == '0:0:0'
     129        assert header['nrofbands'] == '1'
     130        assert header['value'] == '"Default_Band"'
     131
     132        # non-default values
     133        assert header['nroflines'] == '2'
     134        assert header['nrofcellsperline'] == '3'
     135        assert header['registrationcelly'] == '1'
     136
     137
     138    def test_write_non_default_header(self):
     139        header_filename = 'test_write_ermapper_grid.ers'
     140
     141        # setup test data
     142        header = {}
     143        # have some values that aren't defaults
     144        header['nroflines'] = '2'
     145        header['nrofcellsperline'] = '3'
     146       
     147        # Write test header using non-default values                           
     148        ermapper_grids.write_ermapper_header(header_filename, header)
     149
     150        # Check that the read in values match the default values
     151        header = ermapper_grids.read_ermapper_header(header_filename)
     152
     153        # default values
     154        assert header['datum'] == '"GDA94"'
     155        assert header['projection'] == '"GEOGRAPHIC"' 
     156        assert header['coordinatetype'] == 'LL'
     157        assert header['rotation'] == '0:0:0.0'
     158        assert header['celltype'] == 'IEEE4ByteReal'
     159        assert header['nullcellvalue'] == '-99999'
     160        assert header['xdimension'] == '100'
     161        assert header['registrationcellx'] == '0'
     162        assert header['ydimension'] == '100'
     163        assert header['longitude'] == '0:0:0'
     164        assert header['latitude'] == '0:0:0'
     165        assert header['nrofbands'] == '1'
     166        assert header['value'] == '"Default_Band"'
     167
     168        # non-default values
     169        assert header['nroflines'] == '2'
     170        assert header['nrofcellsperline'] == '3'
     171        assert header['registrationcelly'] == '1'
     172
     173        # Clean up created files
     174        remove(header_filename)       
     175
    65176
    66177# def test_default_filenames
Note: See TracChangeset for help on using the changeset viewer.