Python NetCDFFile.variables['elevation'][:]方法代码示例

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

#.........这里部分代码省略.........
    starttime = times[0]

    sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
    sww.store_header(outfile, times, number_of_volumes,
                     number_of_points, description=description,
                     verbose=verbose, sww_precision=netcdf_float)

    # 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')

    # Check zone boundaries
    refzone, _, _ = redfearn(latitudes[0], longitudes[0])

    vertices = {}
    i = 0
    for k, lat in enumerate(latitudes):       # Y direction
        for l, lon in enumerate(longitudes):  # X direction
            vertices[l, k] = i

            _, easting, northing = redfearn(lat, lon)

            #msg = 'Zone boundary crossed at longitude =', lon
            #assert zone == refzone, msg
            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
            x[i] = easting
            y[i] = northing
            i += 1

    #Construct 2 triangles per 'rectangular' element
    volumes = []
    for l in range(number_of_longitudes-1):    # X direction
        for k in range(number_of_latitudes-1): # Y direction
            v1 = vertices[l, k+1]
            v2 = vertices[l, k]
            v3 = vertices[l+1, k+1]
            v4 = vertices[l+1, k]

            volumes.append([v1, v2, v3]) #Upper element
            volumes.append([v4, v3, v2]) #Lower element

    volumes = num.array(volumes, num.int)      #array default#

    if origin is None:
        origin = Geo_reference(refzone, min(x), min(y))
    geo_ref = write_NetCDF_georeference(origin, outfile)

    if elevation is not None:
        z = elevation
    else:
        if inverted_bathymetry:
            z = -1 * elevations
        else:
            z = elevations
    #FIXME: z should be obtained from MOST and passed in here

    #FIXME use the Write_sww instance(sww) to write this info
    z = num.resize(z, outfile.variables['elevation'][:].shape)
    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
    #outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
    outfile.variables['elevation'][:] = z
    outfile.variables['volumes'][:] = volumes.astype(num.int32) #For Opteron 64

    #Time stepping
    stage = outfile.variables['stage']
    xmomentum = outfile.variables['xmomentum']
    ymomentum = outfile.variables['ymomentum']

    if verbose:
        log.critical('Converting quantities')

    n = len(times)
    for j in range(n):
        if verbose and j % ((n+10)/10) == 0:
            log.critical('  Doing %d of %d' % (j, n))

        i = 0
        for k in range(number_of_latitudes):      # Y direction
            for l in range(number_of_longitudes): # X direction
                w = zscale * amplitudes[j, k, l] / 100 + mean_stage
                stage[j, i] = w
                h = w - z[i]
                xmomentum[j, i] = uspeed[j, k, l]/100*h
                ymomentum[j, i] = vspeed[j, k, l]/100*h
                i += 1

    #outfile.close()

    #FIXME: Refactor using code from file_function.statistics
    #Something like print swwstats(swwname)
    if verbose:
        time_info = times, starttime, mint, maxt
        _show_sww_stats(outfile, swwname, geo_ref, time_info)

    outfile.close()

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

#.........这里部分代码省略.........
    number_of_latitudes = latitudes.shape[0]   # Number latitudes
    number_of_longitudes = longitudes.shape[0] # Number longitudes

    # The permutation vector of contains original indices
    # as given in ordering file or None in which case points
    # are assigned the trivial indices enumerating them from
    # 0 to number_of_points-1
    if permutation is None:
        permutation = num.arange(number_of_points, dtype=num.int)

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

    description = 'Converted from URS mux2 files: %s' % basename_in

    # Create new file
    sts = Write_sts()
    sts.store_header(outfile,
                     times+starttime,
                     number_of_points,
                     description=description,
                     verbose=verbose,
                     sts_precision=netcdf_float)

    # 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

    # Check zone boundaries
    if zone is None:
        refzone, _, _ = redfearn(latitudes[0], longitudes[0],
                                 central_meridian=central_meridian)
    else:
        refzone = zone

    old_zone = refzone
    old_easting = 0.0
    old_northing = 0.0

    for i in range(number_of_points):
        computed_zone, easting, northing = redfearn(latitudes[i], longitudes[i],
                                                    zone=zone,
                                                    central_meridian=central_meridian)
        x[i] = easting
        y[i] = northing
        if computed_zone != refzone:
            msg = 'All sts gauges need to be in the same zone. \n'
            msg += 'offending gauge:Zone %d,%.4f, %4f\n' \
                   % (computed_zone, easting, northing)
            msg += 'previous gauge:Zone %d,%.4f, %4f' \
                   % (old_zone, old_easting, old_northing)
            raise Exception, msg
        old_zone = computed_zone
        old_easting = easting
        old_northing = northing

    if origin is None:
        origin = Geo_reference(refzone, min(x), min(y))
    geo_ref = write_NetCDF_georeference(origin, outfile)

    elevation = num.resize(elevation, outfile.variables['elevation'][:].shape)
    outfile.variables['permutation'][:] = permutation.astype(num.int32) # Opteron 64
    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
    outfile.variables['elevation'][:] = elevation

    stage = outfile.variables['stage']
    xmomentum = outfile.variables['xmomentum']
    ymomentum = outfile.variables['ymomentum']

    if verbose: log.critical('Converting quantities')

    for j in range(len(times)):
        for i in range(number_of_points):
            ha = mux['HA'][i,j]
            ua = mux['UA'][i,j]
            va = mux['VA'][i,j]
            if ha == NODATA:
                if verbose:
                    msg = 'Setting nodata value %d to 0 at time = %f, ' \
                          'point = %d' % (ha, times[j], i)
                    log.critical(msg)
                ha = 0.0
                ua = 0.0
                va = 0.0

            w = zscale*ha + mean_stage
            h = w - elevation[i]
            stage[j,i] = w

            xmomentum[j,i] = ua * h
            ymomentum[j,i] = -va * h # South is positive in mux files


    outfile.close()
    
    if verbose:
        log.critical('Wrote sts file ' + stsname)