Python NetCDFFile.createVariable方法代码示例

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import createVariable [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 createVariable [as 别名]
def csv2sts(infile, outfile, latitude = None, longitude = None,
                    verbose = False):
    """
        Take a csv file and convert it to an sts file.
        
        May be used for timeseries, or any other data.
    """
        
    timeseries_data, col_names = load_csv_as_dict(infile, delimiter=' ')
    
    if not col_names:
        raise IOError('csv2sts: file %s is empty or unreadable.' % infile)
    
    if verbose:
        log.critical('csv2sts input data:')
        for col in col_names:
            log.critical('column ' + col + ':')
            log.critical(timeseries_data[col])        

    data_len = len(timeseries_data.values()[0])
    if verbose:
        log.critical('   data length = %d.' % data_len)
    
    fid = NetCDFFile(outfile, netcdf_mode_w)

    fid.createDimension('number_of_timesteps', data_len)

    if latitude:
        fid.latitude = latitude
        
    if longitude:
        fid.longitude = longitude
    
    for col in col_names:
        fid.createVariable(col, netcdf_float, ('number_of_timesteps',))
        
        fid.variables[col][:] = timeseries_data[col]

    fid.close()

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import createVariable [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 createVariable [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 createVariable [as 别名]
def most2nc(input_file, output_file, inverted_bathymetry=False, verbose=True):
    """Convert a MOST file to NetCDF format.

    input_file           the input file to convert
    output_file          the name of the oputput NetCDF file to create
    inverted_bathymetry  ??
    verbose              True if the function is to be verbose
    """

    # variable names
    long_name = 'LON'
    lat_name = 'LAT'
    elev_name = 'ELEVATION'

    # set up bathymetry
    if inverted_bathymetry:
	up = -1.
    else:
	up = +1.

    # read data from the MOST file
    in_file = open(input_file,'r')

    if verbose: log.critical('reading header')

    nx_ny_str = in_file.readline()
    nx_str,ny_str = nx_ny_str.split()
    nx = int(nx_str)
    ny = int(ny_str)
    h1_list=[]
    for i in range(nx):
	h1_list.append(float(in_file.readline()))

    h2_list=[]
    for j in range(ny):
	h2_list.append(float(in_file.readline()))

    h2_list.reverse()

    if verbose: log.critical('reading depths')

    in_depth_list = in_file.readlines()
    in_file.close()

    out_depth_list = [[]]

    if verbose: log.critical('processing depths')

    k=1
    for in_line in in_depth_list:
	for string in in_line.split():
	    #j = k/nx
	    out_depth_list[(k-1)/nx].append(float(string)*up)
	    if k==nx*ny:
		break
	    if k-(k/nx)*nx ==0:
		out_depth_list.append([])
	    k+=1

    in_file.close()
    out_depth_list.reverse()
    depth_list = out_depth_list

    # write the NetCDF file
    if verbose: log.critical('writing results')

    out_file = NetCDFFile(output_file, netcdf_mode_w)

    out_file.createDimension(long_name,nx)

    out_file.createVariable(long_name,'d',(long_name,))
    out_file.variables[long_name].point_spacing='uneven'
    out_file.variables[long_name].units='degrees_east'
    out_file.variables[long_name][:] = h1_list

    out_file.createDimension(lat_name,ny)
    out_file.createVariable(lat_name,'d',(lat_name,))
    out_file.variables[lat_name].point_spacing='uneven'
    out_file.variables[lat_name].units='degrees_north'
    out_file.variables[lat_name][:] = h2_list

    out_file.createVariable(elev_name,'d',(lat_name,long_name))
    out_file.variables[elev_name].point_spacing='uneven'
    out_file.variables[elev_name].units='meters'
    out_file.variables[elev_name][:] = depth_list

    out_file.close()

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

#.........这里部分代码省略.........
        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"))

    # Store
    from anuga.coordinate_transforms.redfearn import redfearn

    x = num.zeros(number_of_points, num.float)  # Easting
    y = num.zeros(number_of_points, num.float)  # Northing

    if verbose:
        log.critical("Making triangular grid")

    # Get zone of 1st point.
    refzone, _, _ = redfearn(latitudes[0], longitudes[0])

    vertices = {}
    i = 0
    for k, lat in enumerate(latitudes):

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

#.........这里部分代码省略.........
            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

        dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
        totalnopoints = nrows * ncols

        # ========================================
        # Do the preceeding with numpy
        # ========================================
        y = num.arange(nrows, dtype=num.float)
        y = yllcorner + (nrows - 1) * cellsize - y * cellsize

        x = num.arange(ncols, dtype=num.float)
        x = xllcorner + x * cellsize

        xx, yy = num.meshgrid(x, y)

        xx = xx.flatten()
        yy = yy.flatten()

        flag = num.logical_and(
            num.logical_and((xx <= easting_max), (xx >= easting_min)),
            num.logical_and((yy <= northing_max), (yy >= northing_min)),
        )

        dem = dem_elevation[:].flatten()

        id = num.where(flag)[0]

        xx = xx[id]
        yy = yy[id]
        dem = dem[id]

        clippednopoints = len(dem)
        # print clippedpoints

        # print xx
        # print yy
        # print dem

        data_flag = dem != NODATA_value

        data_id = num.where(data_flag)

        xx = xx[data_id]
        yy = yy[data_id]
        dem = dem[data_id]

        nn = clippednopoints - len(dem)

        nopoints = len(dem)

        if verbose:
            log.critical("There are %d values in the raster" % totalnopoints)
            log.critical("There are %d values in the clipped raster" % clippednopoints)
            log.critical("There are %d NODATA_values in the clipped raster" % nn)

        outfile.createDimension("number_of_points", nopoints)
        outfile.createDimension("number_of_dimensions", 2)  # This is 2d data

        # Variable definitions
        outfile.createVariable("points", netcdf_float, ("number_of_points", "number_of_dimensions"))
        outfile.createVariable(quantity_name, netcdf_float, ("number_of_points",))

        # Get handles to the variables
        points = outfile.variables["points"]
        elevation = outfile.variables[quantity_name]

        points[:, 0] = xx - easting_min
        points[:, 1] = yy - northing_min
        elevation[:] = dem

        infile.close()
        outfile.close()

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import createVariable [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 createVariable [as 别名]
    def setUp(self):
        import time
        
        self.verbose = Test_File_Conversion.verbose
        # Create basic mesh
        points, vertices, boundary = rectangular(2, 2)

        # Create shallow water domain
        domain = Domain(points, vertices, boundary)
        domain.default_order = 2

        # Set some field values
        domain.set_quantity('elevation', lambda x,y: -x)
        domain.set_quantity('friction', 0.03)


        ######################
        # Boundary conditions
        B = Transmissive_boundary(domain)
        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})


        ######################
        #Initial condition - with jumps
        bed = domain.quantities['elevation'].vertex_values
        stage = num.zeros(bed.shape, num.float)

        h = 0.3
        for i in range(stage.shape[0]):
            if i % 2 == 0:
                stage[i,:] = bed[i,:] + h
            else:
                stage[i,:] = bed[i,:]

        domain.set_quantity('stage', stage)


        domain.distribute_to_vertices_and_edges()               
        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)


        self.domain = domain

        C = domain.get_vertex_coordinates()
        self.X = C[:,0:6:2].copy()
        self.Y = C[:,1:6:2].copy()

        self.F = bed

        #Write A testfile (not realistic. Values aren't realistic)
        self.test_MOST_file = 'most_small'

        longitudes = [150.66667, 150.83334, 151., 151.16667]
        latitudes = [-34.5, -34.33333, -34.16667, -34]

        long_name = 'LON'
        lat_name = 'LAT'

        nx = 4
        ny = 4
        six = 6


        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
            fid = NetCDFFile(self.test_MOST_file + ext, netcdf_mode_w)

            fid.createDimension(long_name,nx)
            fid.createVariable(long_name,netcdf_float,(long_name,))
            fid.variables[long_name].point_spacing='uneven'
            fid.variables[long_name].units='degrees_east'
            fid.variables[long_name][:] = longitudes

            fid.createDimension(lat_name,ny)
            fid.createVariable(lat_name,netcdf_float,(lat_name,))
            fid.variables[lat_name].point_spacing='uneven'
            fid.variables[lat_name].units='degrees_north'
            fid.variables[lat_name][:] = latitudes

            fid.createDimension('TIME',six)
            fid.createVariable('TIME',netcdf_float,('TIME',))
            fid.variables['TIME'].point_spacing='uneven'
            fid.variables['TIME'].units='seconds'
            fid.variables['TIME'][:] = [0.0, 0.1, 0.6, 1.1, 1.6, 2.1]


            name = ext[1:3].upper()
            if name == 'E.': name = 'ELEVATION'
            fid.createVariable(name,netcdf_float,('TIME', lat_name, long_name))
            fid.variables[name].units='CENTIMETERS'
            fid.variables[name].missing_value=-1.e+034

            fid.variables[name][:] = [[[0.3400644, 0, -46.63519, -6.50198],
                                              [-0.1214216, 0, 0, 0],
                                              [0, 0, 0, 0],
                                              [0, 0, 0, 0]],
                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
                                              [0, 2.291054e-005, 2.291054e-005, 0]],
                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
#.........这里部分代码省略.........

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

#.........这里部分代码省略.........
#                #print cellsize
#
#                tpoints[local_index, :] = [x-easting_min, y-northing_min]
#                telev[local_index] = dem_elevation_r[i, j]
#                global_index += 1
#                local_index += 1
#
#        upper_index = global_index
#
#        if upper_index == lower_index + newcols:
#
#            # Seems to be an error with the windows version of
#            # Netcdf. The following gave errors
#
#            try:
#                points[lower_index:upper_index, :] = tpoints
#                elevation[lower_index:upper_index] = telev
#            except:
#                # so used the following if an error occurs
#                for index in range(newcols):
#                    points[index+lower_index, :] = tpoints[index,:]
#                    elevation[index+lower_index] = telev[index]
#
#    assert global_index == nopoints, 'index not equal to number of points'


    #========================================
    # Do the preceeding with numpy
    #========================================
    y = num.arange(nrows,dtype=num.float)
    y = yllcorner + (nrows-1)*cellsize - y*cellsize

    x = num.arange(ncols,dtype=num.float)
    x = xllcorner + x*cellsize

    xx,yy = num.meshgrid(x,y)

    xx = xx.flatten()
    yy = yy.flatten()

    
    flag = num.logical_and(num.logical_and((xx <= easting_max),(xx >= easting_min)),
                           num.logical_and((yy <= northing_max),(yy >= northing_min)))

    
    dem = dem_elevation[:].flatten()


    id = num.where(flag)[0]

    xx = xx[id]
    yy = yy[id]
    dem = dem[id]


    clippednopoints = len(dem)
    #print clippedpoints
    
    #print xx
    #print yy
    #print dem

    data_flag = dem != NODATA_value

    data_id = num.where(data_flag)

    xx = xx[data_id]
    yy = yy[data_id]
    dem = dem[data_id]

    nn = clippednopoints - len(dem)

    nopoints = len(dem)


    if verbose:
        log.critical('There are %d values in the elevation' % totalnopoints)
        log.critical('There are %d values in the clipped elevation'
                     % clippednopoints)
        log.critical('There are %d NODATA_values in the clipped elevation' % nn)

    outfile.createDimension('number_of_points', nopoints)
    outfile.createDimension('number_of_dimensions', 2) #This is 2d data

    # Variable definitions
    outfile.createVariable('points', netcdf_float, ('number_of_points',
                                                    'number_of_dimensions'))
    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))

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

    points[:,0] = xx - easting_min
    points[:,1] = yy - northing_min
    elevation[:] = dem


    infile.close()
    outfile.close()

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import createVariable [as 别名]
    def test_compute_checksum(self):
        """test_compute_checksum(self):

        Check that checksums on files are OK
        """

        from tempfile import mkstemp, mktemp

        # Generate a text file
        tmp_fd , tmp_name = mkstemp(suffix='.tmp', dir='.')
        fid = os.fdopen(tmp_fd, 'w+b')
        string = 'My temp file with textual content. AAAABBBBCCCC1234'
        fid.write(string)
        fid.close()

        # Have to apply the 64 bit fix here since we aren't comparing two
        # files, but rather a string and a file.
        ref_crc = safe_crc(string)

        checksum = compute_checksum(tmp_name)
        assert checksum == ref_crc

        os.remove(tmp_name)



        # Binary file
        tmp_fd , tmp_name = mkstemp(suffix='.tmp', dir='.')
        fid = os.fdopen(tmp_fd, 'w+b')

        string = 'My temp file with binary content. AAAABBBBCCCC1234'
        fid.write(string)
        fid.close()

        ref_crc = safe_crc(string)
        checksum = compute_checksum(tmp_name)

        assert checksum == ref_crc

        os.remove(tmp_name)        

        # Binary NetCDF File X 2 (use mktemp's name)

        try:
            from anuga.file.netcdf import NetCDFFile
        except ImportError:
            # This code is also used by EQRM which does not require NetCDF
            pass
        else:
            test_array = num.array([[7.0, 3.14], [-31.333, 0.0]])

            # First file
            filename1 = mktemp(suffix='.nc', dir='.')
            fid = NetCDFFile(filename1, netcdf_mode_w)
            fid.createDimension('two', 2)
            fid.createVariable('test_array', netcdf_float,
                               ('two', 'two'))
            fid.variables['test_array'][:] = test_array
            fid.close()

            # Second file
            filename2 = mktemp(suffix='.nc', dir='.')
            fid = NetCDFFile(filename2, netcdf_mode_w)
            fid.createDimension('two', 2)
            fid.createVariable('test_array', netcdf_float,
                               ('two', 'two'))
            fid.variables['test_array'][:] = test_array
            fid.close()


            checksum1 = compute_checksum(filename1)
            checksum2 = compute_checksum(filename2)        
            assert checksum1 == checksum2


            os.remove(filename1)
            os.remove(filename2)

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