Python NetCDFFile.description方法代码示例

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import description [as 别名]
def prepare_timeboundary(filename, verbose = False):
    """Convert benchmark 2 time series to NetCDF tms file.
    This is a 'throw-away' code taylor made for files like
    'Benchmark_2_input.txt' from the LWRU2 benchmark
    """

    from anuga.file.netcdf import NetCDFFile

    if verbose: print 'Creating', filename

    # Read the ascii (.txt) version of this file
    fid = open(filename[:-4] + '.txt')

    # Skip first line
    line = fid.readline()

    # Read remaining lines
    lines = fid.readlines()
    fid.close()


    N = len(lines)
    T = num.zeros(N, num.float)  #Time
    Q = num.zeros(N, num.float)  #Values

    for i, line in enumerate(lines):
        fields = line.split()

        T[i] = float(fields[0])
        Q[i] = float(fields[1])


    # Create tms NetCDF file

    fid = NetCDFFile(filename, 'w')
    fid.institution = 'Geoscience Australia'
    fid.description = 'Input wave for Benchmark 2'
    fid.starttime = 0.0
    fid.createDimension('number_of_timesteps', len(T))
    fid.createVariable('time', netcdf_float, ('number_of_timesteps',))
    fid.variables['time'][:] = T

    fid.createVariable('stage', netcdf_float, ('number_of_timesteps',))
    fid.variables['stage'][:] = Q[:]

    fid.createVariable('xmomentum', netcdf_float, ('number_of_timesteps',))
    fid.variables['xmomentum'][:] = 0.0

    fid.createVariable('ymomentum', netcdf_float, ('number_of_timesteps',))
    fid.variables['ymomentum'][:] = 0.0

    fid.close()

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import description [as 别名]
    def helper_write_msh_file(self, filename, l):
        # open the NetCDF file
        fd = NetCDFFile(filename, netcdf_mode_w)
        fd.description = 'Test file - string arrays'

        # convert list of strings to num.array
        al = num.array(string_to_char(l), num.character)

        # write the list
        fd.createDimension('num_of_strings', al.shape[0])
        fd.createDimension('size_of_strings', al.shape[1])

        var = fd.createVariable('strings', netcdf_char,
                                ('num_of_strings', 'size_of_strings'))
        var[:] = al

        fd.close()

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import description [as 别名]

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

    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'

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

    # NetCDF file definition
    fid = NetCDFFile(netcdfname, netcdf_mode_w)

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

    fid.ncols = ncols
    fid.nrows = nrows
    fid.xllcorner = xllcorner
    fid.yllcorner = yllcorner
    fid.cellsize = cellsize
    fid.NODATA_value = NODATA_value

    fid.zone = zone
    fid.false_easting = false_easting
    fid.false_northing = false_northing
    fid.projection = projection
    fid.datum = datum
    fid.units = units

    # dimension definitions
    fid.createDimension('number_of_rows', nrows)
    fid.createDimension('number_of_columns', ncols)

    # variable definitions
    fid.createVariable('elevation', netcdf_float, ('number_of_rows',
                                                   'number_of_columns'))

    # Get handles to the variables
    elevation = fid.variables['elevation']

    #Store data
    import numpy

    datafile = open(name_in)
    elevation[:,:] = numpy.loadtxt(datafile, skiprows=6)
    datafile.close()

#    n = len(lines[6:])
#    for i, line in enumerate(lines[6:]):
#        fields = line.split()
#        if verbose and i % ((n+10)/10) == 0:
#            log.critical('Processing row %d of %d' % (i, nrows))
#
#        if len(fields) != ncols:
#            msg = 'Wrong number of columns in file "%s" line %d\n' % (name_in, i)
#            msg += 'I got %d elements, but there should have been %d\n' % (len(fields), ncols)
#            raise Exception, msg
#
#        elevation[i, :] = num.array([float(x) for x in fields])

    fid.close()

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import description [as 别名]
def _write_msh_file(file_name, mesh):
    """Write .msh NetCDF file

    WARNING: This function mangles the mesh data structure
    """

    # FIXME(Ole and John): We ran into a problem on Bogong (64 bit)
    # where integers appeared as arrays.  This may be similar to
    # problem seen by Steve in changeset:2778 where he had to wrap
    # them in int.  Now we are trying with the native Integer format
    # (Int == 'l' == Int64). However, that caused casting errors, when
    # 64bit arrays are to be assigned to their NetCDF counterparts. It
    # seems that the NetCDF arrays are 32bit even though they are
    # created with the type Int64. Need to look at the NetCDF library
    # in more detail.

    IntType = num.int32
    #IntType = Int

    #print 'mesh vertices',mesh['vertices'].shape


    #the triangulation
    mesh['vertices'] = num.array(mesh['vertices'], num.float)
    mesh['vertex_attribute_titles'] = \
        num.array(string_to_char(mesh['vertex_attribute_titles']), num.character)

    num_attributes = len(mesh['vertex_attribute_titles'])
    num_vertices = mesh['vertices'].shape[0]
    #print 'num_attrib ',num_attributes
    if mesh['vertex_attributes'] != None:
        mesh['vertex_attributes'] = \
            num.array(mesh['vertex_attributes'], num.float)

    if num_attributes > 0 :
        mesh['vertex_attributes'] = \
            num.reshape(mesh['vertex_attributes'],(num_vertices,-1))



    mesh['segments'] = num.array(mesh['segments'], IntType)
    mesh['segment_tags'] = num.array(string_to_char(mesh['segment_tags']),
                                     num.character)
    mesh['triangles'] = num.array(mesh['triangles'], IntType)
    mesh['triangle_tags'] = num.array(string_to_char(mesh['triangle_tags']),
                                      num.character)
    mesh['triangle_neighbors'] = \
        num.array(mesh['triangle_neighbors'], IntType)

    #the outline
    mesh['points'] = num.array(mesh['points'], num.float)
    mesh['point_attributes'] = num.array(mesh['point_attributes'], num.float)
    mesh['outline_segments'] = num.array(mesh['outline_segments'], IntType)
    mesh['outline_segment_tags'] = \
        num.array(string_to_char(mesh['outline_segment_tags']), num.character)
    mesh['holes'] = num.array(mesh['holes'], num.float)
    mesh['regions'] = num.array(mesh['regions'], num.float)
    mesh['region_tags'] = num.array(string_to_char(mesh['region_tags']), num.character)
    mesh['region_max_areas'] = num.array(mesh['region_max_areas'], num.float)

    # NetCDF file definition
    try:
        outfile = NetCDFFile(file_name, netcdf_mode_w)
    except IOError:
        msg = 'File %s could not be created' % file_name
        raise Exception, msg

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

    # dimension definitions - fixed
    outfile.createDimension('num_of_dimensions', 2)     # This is 2d data
    outfile.createDimension('num_of_segment_ends', 2)   # Segs have two points
    outfile.createDimension('num_of_triangle_vertices', 3)
    outfile.createDimension('num_of_triangle_faces', 3)
    outfile.createDimension('num_of_region_max_area', 1)

    # Create dimensions, variables and set the variables

    # trianglulation
    # vertices
    if (mesh['vertices'].shape[0] > 0):
        outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
        outfile.createVariable('vertices', netcdf_float, ('num_of_vertices',
                                                          'num_of_dimensions'))
        outfile.variables['vertices'][:] = mesh['vertices']

        #print 'mesh vertex attributes', mesh['vertex_attributes'].shape
        
        if (mesh['vertex_attributes'] is not None and
            (mesh['vertex_attributes'].shape[0] > 0 and
             mesh['vertex_attributes'].shape[1] > 0)):
            outfile.createDimension('num_of_vertex_attributes',
                                    mesh['vertex_attributes'].shape[1])
            outfile.createDimension('num_of_vertex_attribute_title_chars',
                                    mesh['vertex_attribute_titles'].shape[1])
            outfile.createVariable('vertex_attributes',
                                   netcdf_float,
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import description [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 description [as 别名]

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

    bath_grid = bath_grid[kmin:kmax, lmin:lmax]
    latitudes = latitudes[kmin:kmax]
    longitudes = longitudes[lmin:lmax]
    number_of_latitudes = len(latitudes)
    number_of_longitudes = len(longitudes)
    number_of_times = len(os.listdir(elevation_dir))
    number_of_points = number_of_latitudes * number_of_longitudes
    number_of_volumes = (number_of_latitudes - 1) * (number_of_longitudes - 1) * 2

    # Work out the times
    if len(elevation_files) > 1:
        # Assume: The time period is less than 24hrs.
        time_period = (int(elevation_files[1][-3:]) - int(elevation_files[0][-3:])) * 60 * 60
        times = [x * time_period for x in range(len(elevation_files))]
    else:
        times = [0.0]

    if verbose:
        log.critical("------------------------------------------------")
        log.critical("Statistics:")
        log.critical("  Extent (lat/lon):")
        log.critical("    lat in [%f, %f], len(lat) == %d" % (min(latitudes), max(latitudes), len(latitudes)))
        log.critical("    lon in [%f, %f], len(lon) == %d" % (min(longitudes), max(longitudes), len(longitudes)))
        log.critical("    t in [%f, %f], len(t) == %d" % (min(times), max(times), len(times)))

    ######### WRITE THE SWW FILE #############

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

    # Create new file
    outfile.institution = "Geoscience Australia"
    outfile.description = "Converted from XXX"

    # For sww compatibility
    outfile.smoothing = "Yes"
    outfile.order = 1

    # Start time in seconds since the epoch (midnight 1/1/1970)
    outfile.starttime = starttime = times[0]

    # dimension definitions
    outfile.createDimension("number_of_volumes", number_of_volumes)
    outfile.createDimension("number_of_vertices", 3)
    outfile.createDimension("number_of_points", number_of_points)
    outfile.createDimension("number_of_timesteps", number_of_times)

    # variable definitions
    outfile.createVariable("x", precision, ("number_of_points",))
    outfile.createVariable("y", precision, ("number_of_points",))
    outfile.createVariable("elevation", precision, ("number_of_points",))

    # FIXME: Backwards compatibility
    # outfile.createVariable('z', precision, ('number_of_points',))
    #################################

    outfile.createVariable("volumes", netcdf_int, ("number_of_volumes", "number_of_vertices"))

    outfile.createVariable("time", precision, ("number_of_timesteps",))

    outfile.createVariable("stage", precision, ("number_of_timesteps", "number_of_points"))

    outfile.createVariable("xmomentum", precision, ("number_of_timesteps", "number_of_points"))

    outfile.createVariable("ymomentum", precision, ("number_of_timesteps", "number_of_points"))

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import description [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 description [as 别名]

#.........这里部分代码省略.........
    if file_text[-4:] != '.txt':
        raise IOError('Input file %s should be of type .txt.' % file_text)

    if file_out is None:
        file_out = file_text[:-4] + '.tms'

    fid = open(file_text)
    line = fid.readline()
    fid.close()

    fields = line.split(',')
    msg = "File %s must have the format 'datetime, value0 value1 value2 ...'" \
          % file_text
    assert len(fields) == 2, msg

    if not time_as_seconds:
        try:
            starttime = calendar.timegm(time.strptime(fields[0], time_format))
        except ValueError:
            msg = 'First field in file %s must be' % file_text
            msg += ' date-time with format %s.\n' % time_format
            msg += 'I got %s instead.' % fields[0]
            raise DataTimeError, msg
    else:
        try:
            starttime = float(fields[0])
        except Error:
            msg = "Bad time format"
            raise DataTimeError, msg

    # Split values
    values = []
    for value in fields[1].split():
        values.append(float(value))

    q = ensure_numeric(values)

    msg = 'ERROR: File must contain at least one independent value'
    assert len(q.shape) == 1, msg

    # Read times proper
    from anuga.config import time_format
    import time, calendar

    fid = open(file_text)
    lines = fid.readlines()
    fid.close()

    N = len(lines)
    d = len(q)

    T = num.zeros(N, num.float)       # Time
    Q = num.zeros((N, d), num.float)  # Values

    for i, line in enumerate(lines):
        fields = line.split(',')
        if not time_as_seconds:
            realtime = calendar.timegm(time.strptime(fields[0], time_format))
        else:
            realtime = float(fields[0])
        T[i] = realtime - starttime

        for j, value in enumerate(fields[1].split()):
            Q[i, j] = float(value)

    msg = 'File %s must list time as a monotonuosly ' % file_text
    msg += 'increasing sequence'
    assert num.alltrue(T[1:] - T[:-1] > 0), msg

    #Create NetCDF file
    fid = NetCDFFile(file_out, netcdf_mode_w)

    fid.institution = 'Geoscience Australia'
    fid.description = 'Time series'

    #Reference point
    #Start time in seconds since the epoch (midnight 1/1/1970)
    #FIXME: Use Georef
    fid.starttime = starttime

    # dimension definitions
    #fid.createDimension('number_of_volumes', self.number_of_volumes)
    #fid.createDimension('number_of_vertices', 3)

    fid.createDimension('number_of_timesteps', len(T))

    fid.createVariable('time', netcdf_float, ('number_of_timesteps',))

    fid.variables['time'][:] = T

    for i in range(Q.shape[1]):
        try:
            name = quantity_names[i]
        except:
            name = 'Attribute%d' % i

        fid.createVariable(name, netcdf_float, ('number_of_timesteps',))
        fid.variables[name][:] = Q[:,i]

    fid.close()

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import description [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 description [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")