Changeset 1865

Timestamp:

Oct 4, 2005, 8:39:27 PM (19 years ago)

Author:

ole

Message:

Wrapped sww2asc and sww2ers into one new function: sww2dem
Finished unittest for ers format

Location:

inundation/pyvolution

Files:

• data_manager.py (modified) (5 diffs)
• test_data_manager.py (modified) (8 diffs)

inundation/pyvolution/data_manager.py

@@ -1255 +1255 @@
     false_easting = 500000
     false_northing = 10000000
-    NODATA_value = -9999
 
     if quantity is None:
@@ -1324 +1323 @@
     #Now q has dimension: number_of_points
 
-
+    #Create grid and update xll/yll corner and x,y
+    if verbose: print 'Creating grid'
+    ncols = int((xmax-xmin)/cellsize)+1
+    nrows = int((ymax-ymin)/cellsize)+1
+
+    newxllcorner = xmin+xllcorner
+    newyllcorner = ymin+yllcorner
+
+    x = x+xllcorner-newxllcorner
+    y = y+yllcorner-newyllcorner
+
+    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    assert len(vertex_points.shape) == 2
+
+
+    from Numeric import zeros, Float
+    grid_points = zeros ( (ncols*nrows, 2), Float )
+
+
+    for i in xrange(nrows):
+        yg = i*cellsize
+        for j in xrange(ncols):
+            xg = j*cellsize
+            k = i*ncols + j
+
+            grid_points[k,0] = xg
+            grid_points[k,1] = yg
+
+    #Interpolate
+    from least_squares import Interpolation
+    from util import inside_polygon
+
+    #FIXME: This should be done with precrop = True, otherwise it'll
+    #take forever. With expand_search  set to False, some grid points might
+    #miss out....
+
+    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
+                           precrop = False, expand_search = False,
+                           verbose = verbose)
+
+    #Interpolate using quantity values
+    if verbose: print 'Interpolating'
+    grid_values = interp.interpolate(q).flat
+
+    #Write
     #Write prj file
     if verbose: print 'Writing %s' %prjfile
@@ -1339 +1382 @@
     prjid.close()
 
-    #Create grid and update xll/yll corner and x,y
-    if verbose: print 'Creating grid'
-    ncols = int((xmax-xmin)/cellsize)+1
-    nrows = int((ymax-ymin)/cellsize)+1
-
-    newxllcorner = xmin+xllcorner
-    newyllcorner = ymin+yllcorner
-
-    x = x+xllcorner-newxllcorner
-    y = y+yllcorner-newyllcorner
-
-    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
-    assert len(vertex_points.shape) == 2
-
-
-    from Numeric import zeros, Float
-    grid_points = zeros ( (ncols*nrows, 2), Float )
-
-
-    for i in xrange(nrows):
-        yg = i*cellsize
-        for j in xrange(ncols):
-            xg = j*cellsize
-            k = i*ncols + j
-
-            #print k, xg, yg
-            grid_points[k,0] = xg
-            grid_points[k,1] = yg
-
-
-    #Interpolate
-
-    from least_squares import Interpolation
-    from util import inside_polygon
-
-    #FIXME: This should be done with precrop = True, otherwise it'll
-    #take forever. With expand_search  set to False, some grid points might
-    #miss out....
-
-    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
-                           precrop = False, expand_search = False,
-                           verbose = verbose)
-
-    #Interpolate using quantity values
-    if verbose: print 'Interpolating'
-    grid_values = interp.interpolate(q).flat
-
-    #Write
+
+    
     if verbose: print 'Writing %s' %ascfile
+    NODATA_value = -9999
+ 
     ascid = open(ascfile, 'w')
 
@@ -1419 +1418 @@
     ascid.close()
     fid.close()
+
+def sww2ers(basename_in, basename_out = None,
+            quantity = None,
+            timestep = None,
+            reduction = None,
+            cellsize = 10,
+            verbose = False,
+            origin = None,
+            datum = 'WGS84'):
+    
+    """Read SWW file and convert to Digitial Elevation model format (.asc)
+
+    Example:
+
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+
+    Also write accompanying file with same basename_in but extension .prj
+    used to fix the UTM zone, datum, false northings and eastings.
+
+    The prj format is assumed to be as
+
+    Projection    UTM
+    Zone          56
+    Datum         WGS84
+    Zunits        NO
+    Units         METERS
+    Spheroid      WGS84
+    Xshift        0.0000000000
+    Yshift        10000000.0000000000
+    Parameters
+
+
+    if quantity is given, out values from quantity otherwise default to
+    elevation
+
+    if timestep (an index) is given, output quantity at that timestep
+
+    if reduction is given use that to reduce quantity over all timesteps.
+
+    """
+    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue 
+    import ermapper_grids
+
+    header = {}
+    false_easting = 500000
+    false_northing = 10000000
+    NODATA_value = 0
+
+    if quantity is None:
+        quantity = 'elevation'
+
+    if reduction is None:
+        reduction = max
+
+    if basename_out is None:
+        basename_out = basename_in + '_%s' %quantity
+  
+    swwfile = basename_in + '.sww'
+    # Note the use of a .ers extension is optional (write_ermapper_grid will
+    # deal with either option
+    ersfile = basename_out
+##    prjfile = basename_out + '.prj'
+
+
+    if verbose: print 'Reading from %s' %swwfile
+    #Read sww file
+    from Scientific.IO.NetCDF import NetCDFFile
+    fid = NetCDFFile(swwfile)
+
+    #Get extent and reference
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:]
+
+    ymin = min(y); ymax = max(y)
+    xmin = min(x); xmax = max(x)
+
+    number_of_timesteps = fid.dimensions['number_of_timesteps']
+    number_of_points = fid.dimensions['number_of_points']
+    if origin is None:
+
+        #Get geo_reference
+        #sww files don't have to have a geo_ref
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference() #Default georef object
+
+        xllcorner = geo_reference.get_xllcorner()
+        yllcorner = geo_reference.get_yllcorner()
+        zone = geo_reference.get_zone()
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]
+
+
+    #Get quantity and reduce if applicable
+    if verbose: print 'Reading quantity %s' %quantity
+
+    if quantity.lower() == 'depth':
+        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
+    else:
+        q = fid.variables[quantity][:]
+
+
+    if len(q.shape) == 2:
+        if verbose: print 'Reducing quantity %s' %quantity
+        q_reduced = zeros( number_of_points, Float )
+
+        for k in range(number_of_points):
+            q_reduced[k] = reduction( q[:,k] )
+
+        q = q_reduced
+
+    #Now q has dimension: number_of_points
+
+    #Create grid and update xll/yll corner and x,y
+    if verbose: print 'Creating grid'
+    ncols = int((xmax-xmin)/cellsize)+1
+    nrows = int((ymax-ymin)/cellsize)+1
+
+    newxllcorner = xmin+xllcorner
+    newyllcorner = ymin+yllcorner
+
+    x = x+xllcorner-newxllcorner
+    y = y+yllcorner-newyllcorner
+
+    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    assert len(vertex_points.shape) == 2
+
+
+    from Numeric import zeros, Float
+    grid_points = zeros ( (ncols*nrows, 2), Float )
+
+
+    for i in xrange(nrows):
+        yg = (nrows-i)*cellsize # this will flip the order of the y values
+        for j in xrange(ncols):
+            xg = j*cellsize
+            k = i*ncols + j
+
+            grid_points[k,0] = xg
+            grid_points[k,1] = yg
+
+    #Interpolate
+    from least_squares import Interpolation
+    from util import inside_polygon
+
+    #FIXME: This should be done with precrop = True (?), otherwise it'll
+    #take forever. With expand_search  set to False, some grid points might
+    #miss out....
+
+    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
+                           precrop = False, expand_search = False,
+                           verbose = verbose)
+
+    #Interpolate using quantity values
+    if verbose: print 'Interpolating'
+    grid_values = interp.interpolate(q).flat
+    grid_values = reshape(grid_values,(nrows, ncols))
+
+
+    # setup header information
+    header['datum'] = '"' + datum + '"'
+    # FIXME The use of hardwired UTM and zone number needs to be made optional
+    # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
+    header['projection'] = '"UTM-' + str(zone) + '"'
+    header['coordinatetype'] = 'EN'
+    if header['coordinatetype'] == 'LL':
+        header['longitude'] = str(newxllcorner)
+        header['latitude'] = str(newyllcorner)
+    elif header['coordinatetype'] == 'EN':
+        header['eastings'] = str(newxllcorner)
+        header['northings'] = str(newyllcorner)
+    header['nullcellvalue'] = str(NODATA_value)
+    header['xdimension'] = str(cellsize)
+    header['ydimension'] = str(cellsize)
+    header['value'] = '"' + quantity + '"'
+
+
+    #Write
+    if verbose: print 'Writing %s' %ersfile
+    ermapper_grids.write_ermapper_grid(ersfile,grid_values, header)
+
+    fid.close()    
+
+##    ascid = open(ascfile, 'w')
+##
+##    ascid.write('ncols         %d\n' %ncols)
+##    ascid.write('nrows         %d\n' %nrows)
+##    ascid.write('xllcorner     %d\n' %newxllcorner)
+##    ascid.write('yllcorner     %d\n' %newyllcorner)
+##    ascid.write('cellsize      %f\n' %cellsize)
+##    ascid.write('NODATA_value  %d\n' %NODATA_value)
+##
+##
+##    #Get bounding polygon from mesh
+##    P = interp.mesh.get_boundary_polygon()
+##    inside_indices = inside_polygon(grid_points, P)
+##
+##    for i in range(nrows):
+##        if verbose and i%((nrows+10)/10)==0:
+##            print 'Doing row %d of %d' %(i, nrows)
+##
+##        for j in range(ncols):
+##            index = (nrows-i-1)*ncols+j
+##
+##            if sometrue(inside_indices == index):
+##                ascid.write('%f ' %grid_values[index])
+##            else:
+##                ascid.write('%d ' %NODATA_value)
+##
+##        ascid.write('\n')
+##
+##    #Close
+##    ascid.close()
+##    fid.close()
+
+
+
+def sww2dem(basename_in, basename_out = None,
+            quantity = None,
+            timestep = None,
+            reduction = None,
+            cellsize = 10,
+            verbose = False,
+            origin = None,
+            datum = 'WGS84',
+            format = 'ers'):
+    
+    """Read SWW file and convert to Digitial Elevation model format (.asc or .ers)
+
+    Example (ASC):
+
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+
+    Also write accompanying file with same basename_in but extension .prj
+    used to fix the UTM zone, datum, false northings and eastings.
+
+    The prj format is assumed to be as
+
+    Projection    UTM
+    Zone          56
+    Datum         WGS84
+    Zunits        NO
+    Units         METERS
+    Spheroid      WGS84
+    Xshift        0.0000000000
+    Yshift        10000000.0000000000
+    Parameters
+
+
+    if quantity is given, out values from quantity otherwise default to
+    elevation
+
+    if timestep (an index) is given, output quantity at that timestep
+
+    if reduction is given use that to reduce quantity over all timesteps.
+
+    datum 
+    
+    format can be either 'asc' or 'ers'
+    """
+
+    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue 
+
+    msg = 'Format must be either asc or ers'
+    assert format.lower() in ['asc', 'ers'], msg
+    
+    false_easting = 500000
+    false_northing = 10000000
+
+    if quantity is None:
+        quantity = 'elevation'
+
+    if reduction is None:
+        reduction = max
+
+    if basename_out is None:
+        basename_out = basename_in + '_%s' %quantity
+  
+    swwfile = basename_in + '.sww'
+    demfile = basename_out + '.' + format
+    # Note the use of a .ers extension is optional (write_ermapper_grid will
+    # deal with either option
+    # ersfile = basename_out
+
+    if verbose: print 'Reading from %s' %swwfile
+    #Read sww file
+    from Scientific.IO.NetCDF import NetCDFFile
+    fid = NetCDFFile(swwfile)
+
+    #Get extent and reference
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:]
+
+    ymin = min(y); ymax = max(y)
+    xmin = min(x); xmax = max(x)
+
+    number_of_timesteps = fid.dimensions['number_of_timesteps']
+    number_of_points = fid.dimensions['number_of_points']
+    if origin is None:
+
+        #Get geo_reference
+        #sww files don't have to have a geo_ref
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference() #Default georef object
+
+        xllcorner = geo_reference.get_xllcorner()
+        yllcorner = geo_reference.get_yllcorner()
+        zone = geo_reference.get_zone()
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]
+
+
+    #Get quantity and reduce if applicable
+    if verbose: print 'Reading quantity %s' %quantity
+
+    if quantity.lower() == 'depth':
+        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
+    else:
+        q = fid.variables[quantity][:]
+
+
+    if len(q.shape) == 2:
+        if verbose: print 'Reducing quantity %s' %quantity
+        q_reduced = zeros( number_of_points, Float )
+
+        for k in range(number_of_points):
+            q_reduced[k] = reduction( q[:,k] )
+
+        q = q_reduced
+
+    #Now q has dimension: number_of_points
+
+    #Create grid and update xll/yll corner and x,y
+    if verbose: print 'Creating grid'
+    ncols = int((xmax-xmin)/cellsize)+1
+    nrows = int((ymax-ymin)/cellsize)+1
+
+    newxllcorner = xmin+xllcorner
+    newyllcorner = ymin+yllcorner
+
+    x = x+xllcorner-newxllcorner
+    y = y+yllcorner-newyllcorner
+
+    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    assert len(vertex_points.shape) == 2
+
+
+    from Numeric import zeros, Float
+    grid_points = zeros ( (ncols*nrows, 2), Float )
+
+
+    for i in xrange(nrows):
+        if format.lower() == 'asc':
+            yg = i*cellsize
+        else:    
+            #this will flip the order of the y values for ers
+            yg = (nrows-i)*cellsize 
+            
+        for j in xrange(ncols):
+            xg = j*cellsize
+            k = i*ncols + j
+
+            grid_points[k,0] = xg
+            grid_points[k,1] = yg
+
+    #Interpolate
+    from least_squares import Interpolation
+    from util import inside_polygon
+
+    #FIXME: This should be done with precrop = True (?), otherwise it'll
+    #take forever. With expand_search  set to False, some grid points might
+    #miss out....
+
+    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
+                           precrop = False, expand_search = False,
+                           verbose = verbose)
+
+    #Interpolate using quantity values
+    if verbose: print 'Interpolating'
+    grid_values = interp.interpolate(q).flat
+
+    if format.lower() == 'ers':
+        # setup ERS header information
+        grid_values = reshape(grid_values,(nrows, ncols))    
+        NODATA_value = 0
+        header = {}
+        header['datum'] = '"' + datum + '"'
+        # FIXME The use of hardwired UTM and zone number needs to be made optional
+        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
+        header['projection'] = '"UTM-' + str(zone) + '"'
+        header['coordinatetype'] = 'EN'
+        if header['coordinatetype'] == 'LL':
+            header['longitude'] = str(newxllcorner)
+            header['latitude'] = str(newyllcorner)
+        elif header['coordinatetype'] == 'EN':
+            header['eastings'] = str(newxllcorner)
+            header['northings'] = str(newyllcorner)
+        header['nullcellvalue'] = str(NODATA_value)
+        header['xdimension'] = str(cellsize)
+        header['ydimension'] = str(cellsize)
+        header['value'] = '"' + quantity + '"'
+
+
+        #Write
+        if verbose: print 'Writing %s' %ersfile
+        import ermapper_grids
+        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
+
+        fid.close()    
+    else:
+        #Write to Ascii format
+        
+        #Write prj file
+        prjfile = basename_out + '.prj' 
+        
+        if verbose: print 'Writing %s' %prjfile
+        prjid = open(prjfile, 'w')
+        prjid.write('Projection    %s\n' %'UTM')
+        prjid.write('Zone          %d\n' %zone)
+        prjid.write('Datum         %s\n' %datum)
+        prjid.write('Zunits        NO\n')
+        prjid.write('Units         METERS\n')
+        prjid.write('Spheroid      %s\n' %datum)
+        prjid.write('Xshift        %d\n' %false_easting)
+        prjid.write('Yshift        %d\n' %false_northing)
+        prjid.write('Parameters\n')
+        prjid.close()
+
+
+    
+        if verbose: print 'Writing %s' %ascfile
+        NODATA_value = -9999
+ 
+        ascid = open(demfile, 'w')
+
+        ascid.write('ncols         %d\n' %ncols)
+        ascid.write('nrows         %d\n' %nrows)
+        ascid.write('xllcorner     %d\n' %newxllcorner)
+        ascid.write('yllcorner     %d\n' %newyllcorner)
+        ascid.write('cellsize      %f\n' %cellsize)
+        ascid.write('NODATA_value  %d\n' %NODATA_value)
+
+
+        #Get bounding polygon from mesh
+        P = interp.mesh.get_boundary_polygon()
+        inside_indices = inside_polygon(grid_points, P)
+
+        for i in range(nrows):
+            if verbose and i%((nrows+10)/10)==0:
+                print 'Doing row %d of %d' %(i, nrows)
+
+            for j in range(ncols):
+                index = (nrows-i-1)*ncols+j
+
+                if sometrue(inside_indices == index):
+                    ascid.write('%f ' %grid_values[index])
+                else:
+                    ascid.write('%d ' %NODATA_value)
+
+            ascid.write('\n')
+
+        #Close
+        ascid.close()
+        fid.close()
 
 
@@ -2715 +3199 @@
     outfile.close()
 
-def sww2ers(basename_in, basename_out = None,
-            quantity = None,
-            timestep = None,
-            reduction = None,
-            cellsize = 10,
-            verbose = False,
-            origin = None,
-            datum = 'WGS84'):
-    
-    """Read SWW file and convert to Digitial Elevation model format (.asc)
-
-    Example:
-
-    ncols         3121
-    nrows         1800
-    xllcorner     722000
-    yllcorner     5893000
-    cellsize      25
-    NODATA_value  -9999
-    138.3698 137.4194 136.5062 135.5558 ..........
-
-    Also write accompanying file with same basename_in but extension .prj
-    used to fix the UTM zone, datum, false northings and eastings.
-
-    The prj format is assumed to be as
-
-    Projection    UTM
-    Zone          56
-    Datum         WGS84
-    Zunits        NO
-    Units         METERS
-    Spheroid      WGS84
-    Xshift        0.0000000000
-    Yshift        10000000.0000000000
-    Parameters
-
-
-    if quantity is given, out values from quantity otherwise default to
-    elevation
-
-    if timestep (an index) is given, output quantity at that timestep
-
-    if reduction is given use that to reduce quantity over all timesteps.
-
-    """
-    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue 
-    import ermapper_grids
-
-    header = {}
-    false_easting = 500000
-    false_northing = 10000000
-    NODATA_value = 0
-
-    if quantity is None:
-        quantity = 'elevation'
-
-    if reduction is None:
-        reduction = max
-
-    if basename_out is None:
-        basename_out = basename_in + '_%s' %quantity
-  
-    swwfile = basename_in + '.sww'
-    # Note the use of a .ers extension is optional (write_ermapper_grid will
-    # deal with either option
-    ersfile = basename_out
-##    prjfile = basename_out + '.prj'
-
-
-    if verbose: print 'Reading from %s' %swwfile
-    #Read sww file
-    from Scientific.IO.NetCDF import NetCDFFile
-    fid = NetCDFFile(swwfile)
-
-    #Get extent and reference
-    x = fid.variables['x'][:]
-    y = fid.variables['y'][:]
-    volumes = fid.variables['volumes'][:]
-
-    ymin = min(y); ymax = max(y)
-    xmin = min(x); xmax = max(x)
-
-    number_of_timesteps = fid.dimensions['number_of_timesteps']
-    number_of_points = fid.dimensions['number_of_points']
-    if origin is None:
-
-        #Get geo_reference
-        #sww files don't have to have a geo_ref
-        try:
-            geo_reference = Geo_reference(NetCDFObject=fid)
-        except AttributeError, e:
-            geo_reference = Geo_reference() #Default georef object
-
-        xllcorner = geo_reference.get_xllcorner()
-        yllcorner = geo_reference.get_yllcorner()
-        zone = geo_reference.get_zone()
-    else:
-        zone = origin[0]
-        xllcorner = origin[1]
-        yllcorner = origin[2]
-
-
-    #Get quantity and reduce if applicable
-    if verbose: print 'Reading quantity %s' %quantity
-
-    if quantity.lower() == 'depth':
-        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
-    else:
-        q = fid.variables[quantity][:]
-
-
-    if len(q.shape) == 2:
-        if verbose: print 'Reducing quantity %s' %quantity
-        q_reduced = zeros( number_of_points, Float )
-
-        for k in range(number_of_points):
-            q_reduced[k] = reduction( q[:,k] )
-
-        q = q_reduced
-
-    #Now q has dimension: number_of_points
-
-    #Create grid and update xll/yll corner and x,y
-    if verbose: print 'Creating grid'
-    ncols = int((xmax-xmin)/cellsize)+1
-    nrows = int((ymax-ymin)/cellsize)+1
-
-    newxllcorner = xmin+xllcorner
-    newyllcorner = ymin+yllcorner
-
-    x = x+xllcorner-newxllcorner
-    y = y+yllcorner-newyllcorner
-
-    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
-    assert len(vertex_points.shape) == 2
-
-
-    from Numeric import zeros, Float
-    grid_points = zeros ( (ncols*nrows, 2), Float )
-
-
-    for i in xrange(nrows):
-##        yg = i*cellsize
-        yg = (nrows-i)*cellsize # this will flip the order of the y values
-        for j in xrange(ncols):
-            xg = j*cellsize
-            k = i*ncols + j
-
-##            print k, xg, yg
-            grid_points[k,0] = xg
-            grid_points[k,1] = yg
-
-    #Interpolate
-    from least_squares import Interpolation
-    from util import inside_polygon
-
-    #FIXME: This should be done with precrop = True (?), otherwise it'll
-    #take forever. With expand_search  set to False, some grid points might
-    #miss out....
-
-    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
-                           precrop = False, expand_search = False,
-                           verbose = verbose)
-
-    #Interpolate using quantity values
-    if verbose: print 'Interpolating'
-    grid_values = interp.interpolate(q).flat
-    grid_values = reshape(grid_values,(nrows, ncols))
-
-
-    # setup header information
-    header['datum'] = '"' + datum + '"'
-    # FIXME The use of hardwired UTM and zone number needs to be made optional
-    # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
-    header['projection'] = '"UTM-' + str(zone) + '"'
-    header['coordinatetype'] = 'EN'
-    if header['coordinatetype'] == 'LL':
-        header['longitude'] = str(newxllcorner)
-        header['latitude'] = str(newyllcorner)
-    elif header['coordinatetype'] == 'EN':
-        header['eastings'] = str(newxllcorner)
-        header['northings'] = str(newyllcorner)
-    header['nullcellvalue'] = str(NODATA_value)
-    header['xdimension'] = str(cellsize)
-    header['ydimension'] = str(cellsize)
-    header['value'] = '"' + quantity + '"'
-
-
-    #Write
-    if verbose: print 'Writing %s' %ersfile
-    ermapper_grids.write_ermapper_grid(ersfile,grid_values, header)
-
-    fid.close()    
-
-##    ascid = open(ascfile, 'w')
-##
-##    ascid.write('ncols         %d\n' %ncols)
-##    ascid.write('nrows         %d\n' %nrows)
-##    ascid.write('xllcorner     %d\n' %newxllcorner)
-##    ascid.write('yllcorner     %d\n' %newyllcorner)
-##    ascid.write('cellsize      %f\n' %cellsize)
-##    ascid.write('NODATA_value  %d\n' %NODATA_value)
-##
-##
-##    #Get bounding polygon from mesh
-##    P = interp.mesh.get_boundary_polygon()
-##    inside_indices = inside_polygon(grid_points, P)
-##
-##    for i in range(nrows):
-##        if verbose and i%((nrows+10)/10)==0:
-##            print 'Doing row %d of %d' %(i, nrows)
-##
-##        for j in range(ncols):
-##            index = (nrows-i-1)*ncols+j
-##
-##            if sometrue(inside_indices == index):
-##                ascid.write('%f ' %grid_values[index])
-##            else:
-##                ascid.write('%d ' %NODATA_value)
-##
-##        ascid.write('\n')
-##
-##    #Close
-##    ascid.close()
-##    fid.close()

inundation/pyvolution/test_data_manager.py

@@ -740 +740 @@
 
 
-    def test_sww2asc_elevation(self):
+    def test_sww2dem_asc_elevation(self):
         """Test that sww information can be converted correctly to asc/prj
         format readable by e.g. ArcView
@@ -785 +785 @@
 
         #Export to ascii/prj files
-        sww2asc(self.domain.filename,
+        sww2dem(self.domain.filename,
                 quantity = 'elevation',
                 cellsize = cellsize,
-                verbose = False)
-
-
+                verbose = False,
+                format = 'asc')
+                
         #Check prj (meta data)
         prjid = open(prjfile)
@@ -877 +877 @@
 
 
-    def test_sww2asc_stage_reduction(self):
+    def test_sww2dem_asc_stage_reduction(self):
         """Test that sww information can be converted correctly to asc/prj
         format readable by e.g. ArcView
@@ -925 +925 @@
 
         #Export to ascii/prj files
-        sww2asc(self.domain.filename,
+        sww2dem(self.domain.filename,
                 quantity = 'stage',
                 cellsize = cellsize,
-                reduction = min)
+                reduction = min,
+                format = 'asc')
 
 
@@ -986 +987 @@
 
 
-    def test_sww2asc_missing_points(self):
+    def test_sww2dem_asc_missing_points(self):
         """Test that sww information can be converted correctly to asc/prj
         format readable by e.g. ArcView
@@ -1073 +1074 @@
 
         #Export to ascii/prj files
-        sww2asc(domain.filename,
+        sww2dem(domain.filename,
                 quantity = 'elevation',
                 cellsize = cellsize,
-                verbose = False)
+                verbose = False,
+                format = 'asc')
 
 
@@ -1128 +1130 @@
         os.remove(ascfile)
         os.remove(swwfile)
+
+    def test_sww2ers_simple(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.filename = 'datatest'
+
+        headerfile = self.domain.filename + '.ers'
+        swwfile = self.domain.filename + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep('stage')
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, 'r')
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ers files
+        #sww2ers(self.domain.filename,
+        #        quantity = 'elevation',
+        #        cellsize = cellsize,
+        #        verbose = False)
+                
+        sww2dem(self.domain.filename,
+                quantity = 'elevation',
+                cellsize = cellsize,
+                verbose = False,
+                format = 'ers')
+
+        #Check header data
+        from ermapper_grids import read_ermapper_header, read_ermapper_data
+        
+        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
+        #print header
+        assert header['projection'].lower() == '"utm-56"'
+        assert header['datum'].lower() == '"wgs84"'
+        assert header['units'].lower() == '"meters"'    
+        assert header['value'].lower() == '"elevation"'         
+        assert header['xdimension'] == '0.25'
+        assert header['ydimension'] == '0.25'   
+        assert float(header['eastings']) == 308500.0   #xllcorner
+        assert float(header['northings']) == 6189000.0 #yllcorner       
+        assert int(header['nroflines']) == 5
+        assert int(header['nrofcellsperline']) == 5     
+        assert int(header['nullcellvalue']) == 0   #?           
+        #FIXME - there is more in the header                    
+
+                
+        #Check grid data                
+        grid = read_ermapper_data(self.domain.filename + '_elevation') 
+        
+        #FIXME (Ole): Why is this the desired reference grid for -x-y?
+        ref_grid = [0,      0,     0,     0,     0,
+                    -1,    -1.25, -1.5,  -1.75, -2.0,
+                    -0.75, -1.0,  -1.25, -1.5,  -1.75,              
+                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
+                    -0.25, -0.5,  -0.75, -1.0,  -1.25]              
+                                          
+                                          
+        assert allclose(grid, ref_grid)
+
+        fid.close()
+        
+        #Cleanup
+        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
+        #Done (Ole) - it was because sww2ers didn't close it's sww file 
+        os.remove(sww.filename)
+
+
+    def xxxtestz_sww2ers_real(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+
+
+
+        #Export to ascii/prj files
+        sww2ers('karratha_100m.sww',
+                quantity = 'depth',
+                cellsize = 5,
+                verbose = False)
 
 
@@ -2052 +2167 @@
         os.remove(root + '_100.dem')
 
-    def test_sww2ers_simple(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-        from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile
-
-        #Setup
-        self.domain.filename = 'datatest'
-
-        headerfile = self.domain.filename + '.ers'
-        swwfile = self.domain.filename + '.sww'
-
-        self.domain.set_datadir('.')
-        self.domain.format = 'sww'
-        self.domain.smooth = True
-        self.domain.set_quantity('elevation', lambda x,y: -x-y)
-
-        self.domain.geo_reference = Geo_reference(56,308500,6189000)
-
-        sww = get_dataobject(self.domain)
-        sww.store_connectivity()
-        sww.store_timestep('stage')
-
-        self.domain.evolve_to_end(finaltime = 0.01)
-        sww.store_timestep('stage')
-
-        cellsize = 0.25
-        #Check contents
-        #Get NetCDF
-
-        fid = NetCDFFile(sww.filename, 'r')
-
-        # Get the variables
-        x = fid.variables['x'][:]
-        y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
-        time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:]
-
-
-        #Export to ascii/prj files
-        sww2ers(self.domain.filename,
-                quantity = 'elevation',
-                cellsize = cellsize,
-                verbose = False)
-
-#FIXME (Ole) Why aren't we testing anything
-
-##        #Check prj (meta data)
-##        prjid = open(prjfile)
-##        lines = prjid.readlines()
-##        prjid.close()
-##
-##        L = lines[0].strip().split()
-##        assert L[0].strip().lower() == 'projection'
-##        assert L[1].strip().lower() == 'utm'
-##
-##        L = lines[1].strip().split()
-##        assert L[0].strip().lower() == 'zone'
-##        assert L[1].strip().lower() == '56'
-##
-##        L = lines[2].strip().split()
-##        assert L[0].strip().lower() == 'datum'
-##        assert L[1].strip().lower() == 'wgs84'
-##
-##        L = lines[3].strip().split()
-##        assert L[0].strip().lower() == 'zunits'
-##        assert L[1].strip().lower() == 'no'
-##
-##        L = lines[4].strip().split()
-##        assert L[0].strip().lower() == 'units'
-##        assert L[1].strip().lower() == 'meters'
-##
-##        L = lines[5].strip().split()
-##        assert L[0].strip().lower() == 'spheroid'
-##        assert L[1].strip().lower() == 'wgs84'
-##
-##        L = lines[6].strip().split()
-##        assert L[0].strip().lower() == 'xshift'
-##        assert L[1].strip().lower() == '500000'
-##
-##        L = lines[7].strip().split()
-##        assert L[0].strip().lower() == 'yshift'
-##        assert L[1].strip().lower() == '10000000'
-##
-##        L = lines[8].strip().split()
-##        assert L[0].strip().lower() == 'parameters'
-##
-##
-##        #Check asc file
-##        ascid = open(ascfile)
-##        lines = ascid.readlines()
-##        ascid.close()
-##
-##        L = lines[0].strip().split()
-##        assert L[0].strip().lower() == 'ncols'
-##        assert L[1].strip().lower() == '5'
-##
-##        L = lines[1].strip().split()
-##        assert L[0].strip().lower() == 'nrows'
-##        assert L[1].strip().lower() == '5'
-##
-##        L = lines[2].strip().split()
-##        assert L[0].strip().lower() == 'xllcorner'
-##        assert allclose(float(L[1].strip().lower()), 308500)
-##
-##        L = lines[3].strip().split()
-##        assert L[0].strip().lower() == 'yllcorner'
-##        assert allclose(float(L[1].strip().lower()), 6189000)
-##
-##        L = lines[4].strip().split()
-##        assert L[0].strip().lower() == 'cellsize'
-##        assert allclose(float(L[1].strip().lower()), cellsize)
-##
-##        L = lines[5].strip().split()
-##        assert L[0].strip() == 'NODATA_value'
-##        assert L[1].strip().lower() == '-9999'
-##
-##        #Check grid values
-##        for j in range(5):
-##            L = lines[6+j].strip().split()
-##            y = (4-j) * cellsize
-##            for i in range(5):
-##                assert allclose(float(L[i]), -i*cellsize - y)
-##
-
-        fid.close()
-        
-        #Cleanup
-        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
-        #Done (Ole) - it was because sww2ers didn't close it's sww file 
-        os.remove(sww.filename)
-
-
-    def xxxtestz_sww2ers_real(self):
-        """Test that sww information can be converted correctly to asc/prj
-        format readable by e.g. ArcView
-        """
-
-        import time, os
-        from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile
-
-
-
-
-        #Export to ascii/prj files
-        sww2ers('karratha_100m.sww',
-                quantity = 'depth',
-                cellsize = 5,
-                verbose = False)
-