Python NetCDFFile.yllcorner方法代码示例

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
    def test_read_NetCDFI(self):
        # test if read_NetCDF
        from anuga.file.netcdf import NetCDFFile
        g = Geo_reference(56,1.9,1.9)
        file_name = tempfile.mktemp(".geo_referenceTest")
        
        outfile = NetCDFFile(file_name, netcdf_mode_w)
        outfile.xllcorner = g.get_xllcorner() 
        outfile.yllcorner =  g.get_yllcorner() 
        outfile.zone = g.get_zone()
        outfile.close()
        
        in_file = NetCDFFile(file_name, netcdf_mode_r)
        new_g = Geo_reference(NetCDFObject=in_file)
        in_file.close()
        os.remove(file_name)

        self.assertTrue(g == new_g, ' failed')

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
def _convert_dem_from_ascii2netcdf(name_in, name_out = None,
                                   verbose = False):
    """Read Digital Elevation model from the following ASCII format (.asc)

    Internal function. See public function convert_dem_from_ascii2netcdf
    for details.
    """

    import os
    from anuga.file.netcdf import NetCDFFile

    root = name_in[:-4]

    # Read Meta data
    if verbose: log.critical('Reading METADATA from %s' % (root + '.prj'))

    metadatafile = open(root + '.prj')
    metalines = metadatafile.readlines()
    metadatafile.close()

    L = metalines[0].strip().split()
    assert L[0].strip().lower() == 'projection'
    projection = L[1].strip()                   #TEXT

    L = metalines[1].strip().split()
    assert L[0].strip().lower() == 'zone'
    zone = int(L[1].strip())

    L = metalines[2].strip().split()
    assert L[0].strip().lower() == 'datum'
    datum = L[1].strip()                        #TEXT

    L = metalines[3].strip().split()
    assert L[0].strip().lower() == 'zunits'     #IGNORE
    zunits = L[1].strip()                       #TEXT

    L = metalines[4].strip().split()
    assert L[0].strip().lower() == 'units'
    units = L[1].strip()                        #TEXT

    L = metalines[5].strip().split()
    assert L[0].strip().lower() == 'spheroid'   #IGNORE
    spheroid = L[1].strip()                     #TEXT

    L = metalines[6].strip().split()
    assert L[0].strip().lower() == 'xshift'
    false_easting = float(L[1].strip())

    L = metalines[7].strip().split()
    assert L[0].strip().lower() == 'yshift'
    false_northing = float(L[1].strip())

    if name_in[-4:] != '.asc':
        raise IOError('Input file %s should be of type .asc.' % name_in)

    #Read DEM data
    datafile = open(name_in)

    if verbose: log.critical('Reading DEM from %s' % (name_in))

    lines = datafile.readlines()
    datafile.close()

    if verbose: log.critical('Got %d lines' % len(lines))

    ncols = int(lines[0].split()[1].strip())
    nrows = int(lines[1].split()[1].strip())

    # Do cellsize (line 4) before line 2 and 3
    cellsize = float(lines[4].split()[1].strip())

    # Checks suggested by Joaquim Luis
    # Our internal representation of xllcorner
    # and yllcorner is non-standard.
    xref = lines[2].split()
    if xref[0].strip() == 'xllcorner':
        xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
    elif xref[0].strip() == 'xllcenter':
        xllcorner = float(xref[1].strip())
    else:
        msg = 'Unknown keyword: %s' % xref[0].strip()
        raise Exception, msg

    yref = lines[3].split()
    if yref[0].strip() == 'yllcorner':
        yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
    elif yref[0].strip() == 'yllcenter':
        yllcorner = float(yref[1].strip())
    else:
        msg = 'Unknown keyword: %s' % yref[0].strip()
        raise Exception, msg

    NODATA_value = int(float(lines[5].split()[1].strip()))

    assert len(lines) == nrows + 6

    if name_out == None:
        netcdfname = name_in[:-4]+'.dem'
    else:
        netcdfname = name_out + '.dem'
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
def dem2dem(name_in, stencil, cellsize_new, name_out=None,
                 verbose=False):
    """Read Digitial Elevation model from the following NetCDF format (.dem)

    Example:

    ncols         3121
    nrows         1800
    xllcorner     722000
    yllcorner     5893000
    cellsize      25
    NODATA_value  -9999
    138.3698 137.4194 136.5062 135.5558 ..........

    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
    """

    import os
    from anuga.file.netcdf import NetCDFFile

    if name_in[-4:] != '.dem':
        raise IOError('Input file %s should be of type .dem.' % name_in)

    if name_out != None and basename_out[-4:] != '.dem':
        raise IOError('Input file %s should be of type .dem.' % name_out)

    #Open existing netcdf file to read
    infile = NetCDFFile(name_in, netcdf_mode_r)

    if verbose: log.critical('Reading DEM from %s' % inname)

    # Read metadata (convert from numpy.int32 to int where appropriate)
    ncols = int(infile.ncols)
    nrows = int(infile.nrows)
    xllcorner = infile.xllcorner
    yllcorner = infile.yllcorner
    cellsize = int(infile.cellsize)
    NODATA_value = int(infile.NODATA_value)
    zone = int(infile.zone)
    false_easting = infile.false_easting
    false_northing = infile.false_northing
    projection = infile.projection
    datum = infile.datum
    units = infile.units

    dem_elevation = infile.variables['elevation']

    #Get output file name
    if name_out == None:
        outname = name_in[:-4] + '_' + repr(cellsize_new) + '.dem'
    else:
        outname = name_out

    if verbose: log.critical('Write decimated NetCDF file to %s' % outname)

    #Determine some dimensions for decimated grid
    (nrows_stencil, ncols_stencil) = stencil.shape
    x_offset = ncols_stencil / 2
    y_offset = nrows_stencil / 2
    cellsize_ratio = int(cellsize_new / cellsize)
    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio

    #print type(ncols_new), ncols_new
    
    #Open netcdf file for output
    outfile = NetCDFFile(outname, netcdf_mode_w)

    #Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'NetCDF DEM format for compact and portable ' \
                          'storage of spatial point data'

    #Georeferencing
    outfile.zone = zone
    outfile.projection = projection
    outfile.datum = datum
    outfile.units = units

    outfile.cellsize = cellsize_new
    outfile.NODATA_value = NODATA_value
    outfile.false_easting = false_easting
    outfile.false_northing = false_northing

    outfile.xllcorner = xllcorner + (x_offset * cellsize)
    outfile.yllcorner = yllcorner + (y_offset * cellsize)
    outfile.ncols = ncols_new
    outfile.nrows = nrows_new

    # dimension definition
    #print nrows_new, ncols_new, nrows_new*ncols_new
    #print type(nrows_new), type(ncols_new), type(nrows_new*ncols_new)
    outfile.createDimension('number_of_points', nrows_new*ncols_new)

    # variable definition
    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))

    # Get handle to the variable
    elevation = outfile.variables['elevation']

#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
def _sww_merge_parallel_non_smooth(swwfiles, output,  verbose=False, delete_old=False):
    """
        Merge a list of sww files into a single file.

        Used to merge files created by parallel runs.

        The sww files to be merged must have exactly the same timesteps.

        It is assumed that the separate sww files have been stored in non_smooth
        format.

        Note that some advanced information and custom quantities may not be
        exported.

        swwfiles is a list of .sww files to merge.
        output is the output filename, including .sww extension.
        verbose True to log output information
    """

    if verbose:
        print "MERGING SWW Files"


    first_file = True
    tri_offset = 0
    for filename in swwfiles:
        if verbose:
            print 'Reading file ', filename, ':'

        fid = NetCDFFile(filename, netcdf_mode_r)

        if first_file:

            times    = fid.variables['time'][:]
            n_steps = len(times)
            number_of_timesteps = fid.dimensions['number_of_timesteps']
            #print n_steps, number_of_timesteps
            starttime = int(fid.starttime)

            out_s_quantities = {}
            out_d_quantities = {}

            out_s_c_quantities = {}
            out_d_c_quantities = {}


            xllcorner = fid.xllcorner
            yllcorner = fid.yllcorner

            number_of_global_triangles = int(fid.number_of_global_triangles)
            number_of_global_nodes     = int(fid.number_of_global_nodes)
            number_of_global_triangle_vertices = 3*number_of_global_triangles


            order      = fid.order
            xllcorner  = fid.xllcorner;
            yllcorner  = fid.yllcorner ;
            zone       = fid.zone;
            false_easting  = fid.false_easting;
            false_northing = fid.false_northing;
            datum      = fid.datum;
            projection = fid.projection;

            g_volumes = num.arange(number_of_global_triangles*3).reshape(-1,3)



            g_x = num.zeros((number_of_global_triangle_vertices,),num.float32)
            g_y = num.zeros((number_of_global_triangle_vertices,),num.float32)

            g_points = num.zeros((number_of_global_triangle_vertices,2),num.float32)

            #=======================================
            # Deal with the vertex based variables
            #=======================================
            quantities = set(['elevation', 'friction', 'stage', 'xmomentum',
                              'ymomentum', 'xvelocity', 'yvelocity', 'height'])
            variables = set(fid.variables.keys())

            quantities = list(quantities & variables)

            static_quantities = []
            dynamic_quantities = []

            for quantity in quantities:
                # Test if elevation is static
                if n_steps == fid.variables[quantity].shape[0]:
                    dynamic_quantities.append(quantity)
                else:
                    static_quantities.append(quantity)

            # Static Quantities are stored as a 1D array
            for quantity in static_quantities:
                out_s_quantities[quantity] = num.zeros((3*number_of_global_triangles,),num.float32)

            #=======================================
            # Deal with the centroid based variables
            #=======================================
            quantities = set(['elevation_c', 'friction_c', 'stage_c', 'xmomentum_c',
                              'ymomentum_c', 'xvelocity_c', 'yvelocity_c', 'height_c'])
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
def _sww_merge(swwfiles, output, verbose=False):
    """
        Merge a list of sww files into a single file.
        
        May be useful for parallel runs. Note that colinear points and
        edges are not merged: there will essentially be multiple meshes within
        the one sww file.
        
        The sww files to be merged must have exactly the same timesteps. Note
        that some advanced information and custom quantities may not be
        exported.
        
        swwfiles is a list of .sww files to merge.
        output is the output filename, including .sww extension.
        verbose True to log output information
    """

    if verbose:
        print "MERGING SWW Files"
        
    static_quantities = ['elevation']
    dynamic_quantities = ['stage', 'xmomentum', 'ymomentum']
    
    first_file = True
    tri_offset = 0
    for filename in swwfiles:
        if verbose:
            print 'Reading file ', filename, ':'    
    
        fid = NetCDFFile(filename, netcdf_mode_r)



        
        
        tris = fid.variables['volumes'][:]       
         
        if first_file:
            times = fid.variables['time'][:]
            x = []
            y = []
            out_tris = list(tris)  
            out_s_quantities = {}
            out_d_quantities = {}


            xllcorner = fid.xllcorner
            yllcorner = fid.yllcorner

            order      = fid.order
            xllcorner  = fid.xllcorner;
            yllcorner  = fid.yllcorner ;
            zone       = fid.zone;
            false_easting  = fid.false_easting;
            false_northing = fid.false_northing;
            datum      = fid.datum;
            projection = fid.projection;

            
            for quantity in static_quantities:
                out_s_quantities[quantity] = []

            # Quantities are stored as a 2D array of timesteps x data.
            for quantity in dynamic_quantities:
                out_d_quantities[quantity] = [ [] for _ in range(len(times))]
                 
            description = 'merged:' + getattr(fid, 'description')          
            first_file = False
        else:
            for tri in tris:
                # Advance new tri indices to point at newly appended points.
                verts = [vertex+tri_offset for vertex in tri]
                out_tris.append(verts)



        try: # works with netcdf4
            num_pts = len(fid.dimensions['number_of_points'])
        except: # works with scientific.io.netcdf
            num_pts = int(fid.dimensions['number_of_points'])

        tri_offset += num_pts
        
        if verbose:
            print '  new triangle index offset is ', tri_offset
            
        x.extend(list(fid.variables['x'][:]))
        y.extend(list(fid.variables['y'][:]))
        
        # Grow the list of static quantities associated with the x,y points
        for quantity in static_quantities:
            out_s_quantities[quantity].extend(fid.variables[quantity][:])
            
        #Collate all dynamic quantities according to their timestep
        for quantity in dynamic_quantities:
            time_chunks = fid.variables[quantity][:]
            for i, time_chunk in enumerate(time_chunks):
                out_d_quantities[quantity][i].extend(time_chunk)            
    
    # Mash all points into a single big list    
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
def _sww_merge_parallel_smooth(swwfiles, output,  verbose=False, delete_old=False):
    """
        Merge a list of sww files into a single file.
        
        Use to merge files created by parallel runs.

        The sww files to be merged must have exactly the same timesteps.

        It is assumed that the separate sww files have been stored in non_smooth
        format.

        Note that some advanced information and custom quantities may not be
        exported.
        
        swwfiles is a list of .sww files to merge.
        output is the output filename, including .sww extension.
        verbose True to log output information
    """

    if verbose:
        print "MERGING SWW Files"
        
    
    first_file = True
    tri_offset = 0
    for filename in swwfiles:
        if verbose:
            print 'Reading file ', filename, ':'    
    
        fid = NetCDFFile(filename, netcdf_mode_r)
         
        if first_file:

            times    = fid.variables['time'][:]
            n_steps = len(times)
            #number_of_timesteps = fid.dimensions['number_of_timesteps']
            #print n_steps, number_of_timesteps
            starttime = int(fid.starttime)
            
            out_s_quantities = {}
            out_d_quantities = {}

            out_s_c_quantities = {}
            out_d_c_quantities = {}


            xllcorner = fid.xllcorner
            yllcorner = fid.yllcorner

            number_of_global_triangles = int(fid.number_of_global_triangles)
            number_of_global_nodes     = int(fid.number_of_global_nodes)

            order      = fid.order
            xllcorner  = fid.xllcorner;
            yllcorner  = fid.yllcorner ;
            zone       = fid.zone;
            false_easting  = fid.false_easting;
            false_northing = fid.false_northing;
            datum      = fid.datum;
            projection = fid.projection;

            g_volumes = num.zeros((number_of_global_triangles,3),num.int)
            g_x = num.zeros((number_of_global_nodes,),num.float32)
            g_y = num.zeros((number_of_global_nodes,),num.float32)

            g_points = num.zeros((number_of_global_nodes,2),num.float32)

            #=====================================
            # Deal with the vertex based variables
            #=====================================
            quantities = set(['elevation', 'friction', 'stage', 'xmomentum',
                              'ymomentum', 'xvelocity', 'yvelocity', 'height'])
            variables = set(fid.variables.keys())

            quantities = list(quantities & variables)
            
            static_quantities = []
            dynamic_quantities = []

            for quantity in quantities:
                # Test if quantity is static
                if n_steps == fid.variables[quantity].shape[0]:
                    dynamic_quantities.append(quantity)
                else:
                    static_quantities.append(quantity)
                
            for quantity in static_quantities:
                out_s_quantities[quantity] = num.zeros((number_of_global_nodes,),num.float32)

            # Quantities are stored as a 2D array of timesteps x data.
            for quantity in dynamic_quantities:
                out_d_quantities[quantity] = \
                      num.zeros((n_steps,number_of_global_nodes),num.float32)

            #=======================================
            # Deal with the centroid based variables
            #=======================================
            quantities = set(['elevation_c', 'friction_c', 'stage_c', 'xmomentum_c',
                              'ymomentum_c', 'xvelocity_c', 'yvelocity_c', 'height_c'])
            variables = set(fid.variables.keys())
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
    def _generic_dem2pts(
        self,
        name_in,
        name_out=None,
        quantity_name=None,
        verbose=False,
        easting_min=None,
        easting_max=None,
        northing_min=None,
        northing_max=None,
    ):
        """Read raster from the following NetCDF format (.dem)

        Internal function. See public function generic_dem2pts for details.
        """

        # FIXME: Can this be written feasibly using write_pts?

        import os
        from anuga.file.netcdf import NetCDFFile

        root = name_in[:-4]

        if name_in[-4:] == ".asc":
            intermediate = root + ".dem"
            if verbose:
                log.critical("Preconvert %s from asc to %s" % (name_in, intermediate))
            asc2dem(name_in)
            name_in = intermediate
        elif name_in[-4:] != ".dem":
            raise IOError("Input file %s should be of type .asc or .dem." % name_in)

        if name_out != None and basename_out[-4:] != ".pts":
            raise IOError("Input file %s should be of type .pts." % name_out)

        # Get NetCDF
        infile = NetCDFFile(name_in, netcdf_mode_r)

        if verbose:
            log.critical("Reading raster from %s" % (name_in))

        ncols = int(infile.ncols)
        nrows = int(infile.nrows)
        xllcorner = float(infile.xllcorner)  # Easting of lower left corner
        yllcorner = float(infile.yllcorner)  # Northing of lower left corner
        cellsize = float(infile.cellsize)
        NODATA_value = float(infile.NODATA_value)

        dem_elevation = infile.variables[quantity_name]

        zone = int(infile.zone)
        false_easting = float(infile.false_easting)
        false_northing = float(infile.false_northing)

        # print ncols, nrows, xllcorner,yllcorner, cellsize, NODATA_value, zone

        # Text strings
        projection = infile.projection
        datum = infile.datum
        units = infile.units

        # print projection, datum, units

        # Get output file
        if name_out == None:
            ptsname = root + ".pts"
        else:
            ptsname = name_out

        if verbose:
            log.critical("Store to NetCDF file %s" % ptsname)

        # NetCDF file definition
        outfile = NetCDFFile(ptsname, netcdf_mode_w)

        # Create new file
        outfile.institution = "Geoscience Australia"
        outfile.description = "NetCDF pts format for compact and portable " "storage of spatial point data"

        # Assign default values
        if easting_min is None:
            easting_min = xllcorner
        if easting_max is None:
            easting_max = xllcorner + ncols * cellsize
        if northing_min is None:
            northing_min = yllcorner
        if northing_max is None:
            northing_max = yllcorner + nrows * cellsize

        # print easting_min, easting_max, northing_min, northing_max

        # Compute offsets to update georeferencing
        easting_offset = xllcorner - easting_min
        northing_offset = yllcorner - northing_min

        # Georeferencing
        outfile.zone = zone
        outfile.xllcorner = easting_min  # Easting of lower left corner
        outfile.yllcorner = northing_min  # Northing of lower left corner
        outfile.false_easting = false_easting
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
def _dem2pts(name_in, name_out=None, verbose=False,
            easting_min=None, easting_max=None,
            northing_min=None, northing_max=None):
    """Read Digitial Elevation model from the following NetCDF format (.dem)

    Internal function. See public function dem2pts for details.
    """

    # FIXME: Can this be written feasibly using write_pts?

    import os
    from anuga.file.netcdf import NetCDFFile

    root = name_in[:-4]

    if name_in[-4:] == '.asc':
        intermediate = root + '.dem'
        if verbose:
            log.critical('Preconvert %s from asc to %s' % \
                                    (name_in, intermediate))
        asc2dem(name_in)
        name_in = intermediate
    elif name_in[-4:] != '.dem':
        raise IOError('Input file %s should be of type .asc or .dem.' % name_in)

    if name_out != None and basename_out[-4:] != '.pts':
        raise IOError('Input file %s should be of type .pts.' % name_out)

    # Get NetCDF
    infile = NetCDFFile(name_in, netcdf_mode_r) 

    if verbose: log.critical('Reading DEM from %s' % (name_in))

    ncols = int(infile.ncols)
    nrows = int(infile.nrows)
    xllcorner = float(infile.xllcorner)  # Easting of lower left corner
    yllcorner = float(infile.yllcorner)  # Northing of lower left corner
    cellsize = float(infile.cellsize)
    NODATA_value = float(infile.NODATA_value)

    dem_elevation = infile.variables['elevation']

    zone = int(infile.zone)
    false_easting = float(infile.false_easting)
    false_northing = float(infile.false_northing)

    #print ncols, nrows, xllcorner,yllcorner, cellsize, NODATA_value, zone


    # Text strings
    projection = infile.projection
    datum = infile.datum
    units = infile.units

    #print projection, datum, units

    # Get output file
    if name_out == None:
        ptsname = root + '.pts'
    else:
        ptsname = name_out

    if verbose: log.critical('Store to NetCDF file %s' % ptsname)

    # NetCDF file definition
    outfile = NetCDFFile(ptsname, netcdf_mode_w)

    # Create new file
    outfile.institution = 'Geoscience Australia'
    outfile.description = 'NetCDF pts format for compact and portable ' \
                          'storage of spatial point data'

    # Assign default values
    if easting_min is None: easting_min = xllcorner
    if easting_max is None: easting_max = xllcorner + ncols*cellsize
    if northing_min is None: northing_min = yllcorner
    if northing_max is None: northing_max = yllcorner + nrows*cellsize


    #print easting_min, easting_max, northing_min, northing_max

    # Compute offsets to update georeferencing
    easting_offset = xllcorner - easting_min
    northing_offset = yllcorner - northing_min

    # Georeferencing
    outfile.zone = zone
    outfile.xllcorner = easting_min # Easting of lower left corner
    outfile.yllcorner = northing_min # Northing of lower left corner
    outfile.false_easting = false_easting
    outfile.false_northing = false_northing

    outfile.projection = projection
    outfile.datum = datum
    outfile.units = units

    # Grid info (FIXME: probably not going to be used, but heck)
    outfile.ncols = ncols
    outfile.nrows = nrows

#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import yllcorner [as 别名]
    def test_decimate_dem(self):
        """Test decimation of dem file
        """

        import os
        from anuga.file.netcdf import NetCDFFile

        # Write test dem file
        root = "decdemtest"

        filename = root + ".dem"
        fid = NetCDFFile(filename, netcdf_mode_w)

        fid.institution = "Geoscience Australia"
        fid.description = "NetCDF DEM format for compact and portable " + "storage of spatial point data"

        nrows = 15
        ncols = 18

        fid.ncols = ncols
        fid.nrows = nrows
        fid.xllcorner = 2000.5
        fid.yllcorner = 3000.5
        fid.cellsize = 25
        fid.NODATA_value = -9999

        fid.zone = 56
        fid.false_easting = 0.0
        fid.false_northing = 0.0
        fid.projection = "UTM"
        fid.datum = "WGS84"
        fid.units = "METERS"

        fid.createDimension("number_of_points", nrows * ncols)

        fid.createVariable("elevation", netcdf_float, ("number_of_points",))

        elevation = fid.variables["elevation"]

        elevation[:] = num.arange(nrows * ncols)

        fid.close()

        # generate the elevation values expected in the decimated file
        ref_elevation = [
            (0 + 1 + 2 + 18 + 19 + 20 + 36 + 37 + 38) / 9.0,
            (4 + 5 + 6 + 22 + 23 + 24 + 40 + 41 + 42) / 9.0,
            (8 + 9 + 10 + 26 + 27 + 28 + 44 + 45 + 46) / 9.0,
            (12 + 13 + 14 + 30 + 31 + 32 + 48 + 49 + 50) / 9.0,
            (72 + 73 + 74 + 90 + 91 + 92 + 108 + 109 + 110) / 9.0,
            (76 + 77 + 78 + 94 + 95 + 96 + 112 + 113 + 114) / 9.0,
            (80 + 81 + 82 + 98 + 99 + 100 + 116 + 117 + 118) / 9.0,
            (84 + 85 + 86 + 102 + 103 + 104 + 120 + 121 + 122) / 9.0,
            (144 + 145 + 146 + 162 + 163 + 164 + 180 + 181 + 182) / 9.0,
            (148 + 149 + 150 + 166 + 167 + 168 + 184 + 185 + 186) / 9.0,
            (152 + 153 + 154 + 170 + 171 + 172 + 188 + 189 + 190) / 9.0,
            (156 + 157 + 158 + 174 + 175 + 176 + 192 + 193 + 194) / 9.0,
            (216 + 217 + 218 + 234 + 235 + 236 + 252 + 253 + 254) / 9.0,
            (220 + 221 + 222 + 238 + 239 + 240 + 256 + 257 + 258) / 9.0,
            (224 + 225 + 226 + 242 + 243 + 244 + 260 + 261 + 262) / 9.0,
            (228 + 229 + 230 + 246 + 247 + 248 + 264 + 265 + 266) / 9.0,
        ]

        # generate a stencil for computing the decimated values
        stencil = num.ones((3, 3), num.float) / 9.0

        dem2dem(filename, stencil=stencil, cellsize_new=100)

        # Open decimated NetCDF file
        fid = NetCDFFile(root + "_100.dem", netcdf_mode_r)

        # Get decimated elevation
        elevation = fid.variables["elevation"]

        # Check values
        assert num.allclose(elevation, ref_elevation)

        # Cleanup
        fid.close()

        os.remove(root + ".dem")
        os.remove(root + "_100.dem")