Python NetCDFFile.zscale方法代码示例

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import zscale [as 别名]
def sts2sww_mesh(basename_in, basename_out=None, 
                 spatial_thinning=1, verbose=False):
    
    from anuga.mesh_engine.mesh_engine import NoTrianglesError
    from anuga.pmesh.mesh import Mesh
    if verbose:
        print "Starting sts2sww_mesh"
    
    mean_stage=0.
    zscale=1.

    if (basename_in[:-4]=='.sts'):
        stsname = basename_in
    else: 
        stsname = basename_in + '.sts'

    if verbose: print "Reading sts NetCDF file: %s" %stsname
    infile = NetCDFFile(stsname, netcdf_mode_r)
    cellsize = infile.cellsize
    ncols = infile.ncols
    nrows = infile.nrows
    no_data = infile.no_data
    refzone = infile.zone
    x_origin = infile.xllcorner
    y_origin = infile.yllcorner
    origin = num.array([x_origin, y_origin])
    x = infile.variables['x'][:]
    y = infile.variables['y'][:]
    times = infile.variables['time'][:]
    wind_speed_full = infile.variables['wind_speed'][:]
    wind_angle_full = infile.variables['wind_angle'][:]
    pressure_full   =   infile.variables['barometric_pressure'][:]
    infile.close()

    number_of_points = nrows*ncols
    points_utm = num.zeros((number_of_points,2),num.float)
    points_utm[:,0]=x+x_origin
    points_utm[:,1]=y+y_origin

    thinned_indices=[]
    for i in range(number_of_points):
        if (i/ncols==0 or i/ncols==ncols-1 or (i/ncols)%(spatial_thinning)==0):
            if ( i%(spatial_thinning)==0 or i%nrows==0 or i%nrows==nrows-1 ):  
                thinned_indices.append(i)

    #Spatial thinning
    points_utm=points_utm[thinned_indices]
    number_of_points = points_utm.shape[0]
    number_of_timesteps = wind_speed_full.shape[0]
    wind_speed = num.empty((number_of_timesteps,number_of_points),dtype=float)
    wind_angle = num.empty((number_of_timesteps,number_of_points),dtype=float)
    barometric_pressure   = num.empty((number_of_timesteps,number_of_points),dtype=float)
    if verbose:
        print "Total number of points: ", nrows*ncols
        print "Number of thinned points: ", number_of_points
    for i in xrange(number_of_timesteps):
        wind_speed[i] = wind_speed_full[i,thinned_indices]
        wind_angle[i] = wind_angle_full[i,thinned_indices]
        barometric_pressure[i]   = pressure_full[i,thinned_indices]

    #P.plot(points_utm[:,0],points_utm[:,1],'ro')
    #P.show()

    if verbose:
        print "Generating sww triangulation of gems data"

    mesh = Mesh()
    mesh.add_vertices(points_utm)
    mesh.auto_segment(smooth_indents=True, expand_pinch=True)
    mesh.auto_segment(mesh.shape.get_alpha() * 1.1)
    try:
        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)
    except NoTrianglesError:
        # This is a bit of a hack, going in and changing the data structure.
        mesh.holes = []
        mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False)

    mesh_dic = mesh.Mesh2MeshList()

    points_utm=ensure_numeric(points_utm)
    assert num.alltrue(ensure_numeric(mesh_dic['generatedpointlist'])
                       == ensure_numeric(points_utm))

    volumes = mesh_dic['generatedtrianglelist']

    # Write sww intro and grid stuff.
    if (basename_out is not None and basename_out[:-4]=='.sww'): 
        swwname = basename_out
    else: 
        swwname = basename_in + '.sww'

    if verbose: 'Output to %s' % swwname

    if verbose:
        print "Writing sww wind and pressure field file"
    outfile = NetCDFFile(swwname, netcdf_mode_w)
    sww = Write_sww([], ['wind_speed','wind_angle','barometric_pressure'])
    sww.store_header(outfile, times, len(volumes), len(points_utm),
                     verbose=verbose, sww_precision='d')
    outfile.mean_stage = mean_stage
#.........这里部分代码省略.........

# 需要导入模块: from anuga.file.netcdf import NetCDFFile [as 别名]
# 或者: from anuga.file.netcdf.NetCDFFile import zscale [as 别名]
def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
                      mint=None, maxt=None,
                      mean_stage=0,
                      origin=None,
                      hole_points_UTM=None,
                      zscale=1):
    """
    Convert URS C binary format for wave propagation to
    sww format native to abstract_2d_finite_volumes.

    Specify only basename_in and read files of the form
    basefilename-z-mux, basefilename-e-mux and
    basefilename-n-mux containing relative height,
    x-velocity and y-velocity, respectively.

    Also convert latitude and longitude to UTM. All coordinates are
    assumed to be given in the GDA94 datum. The latitude and longitude
    information is assumed ungridded grid.

    min's and max's: If omitted - full extend is used.
    To include a value min ans max may equal it.
    Lat and lon are assumed to be in decimal degrees.

    origin is a 3-tuple with geo referenced
    UTM coordinates (zone, easting, northing)
    It will be the origin of the sww file. This shouldn't be used,
    since all of anuga should be able to handle an arbitary origin.
    The mux point info is NOT relative to this origin.

    URS C binary format has data organised as TIME, LONGITUDE, LATITUDE
    which means that latitude is the fastest
    varying dimension (row major order, so to speak)

    In URS C binary the latitudes and longitudes are in assending order.

    Note, interpolations of the resulting sww file will be different
    from results of urs2sww.  This is due to the interpolation
    function used, and the different grid structure between urs2sww
    and this function.

    Interpolating data that has an underlying gridded source can
    easily end up with different values, depending on the underlying
    mesh.

    consider these 4 points
    50  -50

    0     0

    The grid can be
     -
    |\|   A
     -
     or;
      -
     |/|  B
      -

    If a point is just below the center of the midpoint, it will have a
    +ve value in grid A and a -ve value in grid B.
    """

    from anuga.mesh_engine.mesh_engine import NoTrianglesError
    from anuga.pmesh.mesh import Mesh

    files_in = [basename_in + WAVEHEIGHT_MUX_LABEL,
                basename_in + EAST_VELOCITY_LABEL,
                basename_in + NORTH_VELOCITY_LABEL]
    quantities = ['HA','UA','VA']

    # instantiate urs_points of the three mux files.
    mux = {}
    for quantity, file in map(None, quantities, files_in):
        mux[quantity] = Read_urs(file)

    # Could check that the depth is the same. (hashing)

    # handle to a mux file to do depth stuff
    a_mux = mux[quantities[0]]

    # Convert to utm
    lat = a_mux.lonlatdep[:,1]
    long = a_mux.lonlatdep[:,0]
    points_utm, zone = convert_from_latlon_to_utm(latitudes=lat,
                                                  longitudes=long)

    elevation = a_mux.lonlatdep[:,2] * -1

    # grid (create a mesh from the selected points)
    # This mesh has a problem.  Triangles are streched over ungridded areas.
    # If these areas could be described as holes in pmesh, that would be great.

    # I can't just get the user to selection a point in the middle.
    # A boundary is needed around these points.
    # But if the zone of points is obvious enough auto-segment should do
    # a good boundary.
    mesh = Mesh()
    mesh.add_vertices(points_utm)
    mesh.auto_segment(smooth_indents=True, expand_pinch=True)

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