Python NetCDFFile.variables["x"][:]方法代码示例

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

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

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

            zone, 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]

            # Note, this is different to the ferrit2sww code
            # since the order of the lats is reversed.
            volumes.append([v1, v3, v2])  # Upper element
            volumes.append([v4, v2, v3])  # Lower element

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

    geo_ref = Geo_reference(refzone, min(x), min(y))
    geo_ref.write_NetCDF(outfile)

    # This will put the geo ref in the middle
    # geo_ref = Geo_reference(refzone, (max(x)+min(x))/2., (max(x)+min(y))/2.)

    if verbose:
        log.critical("------------------------------------------------")
        log.critical("More Statistics:")
        log.critical("  Extent (/lon):")
        log.critical("    x in [%f, %f], len(lat) == %d" % (min(x), max(x), len(x)))
        log.critical("    y in [%f, %f], len(lon) == %d" % (min(y), max(y), len(y)))
        log.critical("geo_ref: ", geo_ref)

    z = num.resize(bath_grid, outfile.variables["elevation"][:].shape)
    outfile.variables["x"][:] = x - geo_ref.get_xllcorner()
    outfile.variables["y"][:] = y - geo_ref.get_yllcorner()
    # FIXME (Ole): Remove once viewer has been recompiled and changed
    #              to use elevation instead of z
    # outfile.variables['z'][:] = z
    outfile.variables["elevation"][:] = z
    outfile.variables["volumes"][:] = volumes.astype(num.int32)  # On Opteron 64

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

    outfile.variables["time"][:] = times  # Store time relative

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

    n = number_of_times
    for j in range(number_of_times):
        # load in files
        elevation_meta, elevation_grid = _read_asc(elevation_dir + os.sep + elevation_files[j])

        _, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j])
        _, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j])

        # cut matrix to desired size
        elevation_grid = elevation_grid[kmin:kmax, lmin:lmax]
        u_momentum_grid = u_momentum_grid[kmin:kmax, lmin:lmax]
        v_momentum_grid = v_momentum_grid[kmin:kmax, lmin:lmax]

        # handle missing values
        missing = elevation_grid == elevation_meta["NODATA_value"]
        if num.sometrue(missing):
            if fail_on_NaN:
                msg = "File %s contains missing values" % (elevation_files[j])
                raise DataMissingValuesError, msg
            else:
                elevation_grid = elevation_grid * (missing == 0) + missing * elevation_NaN_filler

        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 * elevation_grid[k, l] + mean_stage
                stage[j, i] = w
                h = w - z[i]
                xmomentum[j, i] = u_momentum_grid[k, l] * h
                ymomentum[j, i] = v_momentum_grid[k, l] * h
                i += 1

    outfile.close()