Trond,
See if the attached file contains something close to what you need. It
has no loops at all; I have not timed it, but it should be quite quick.
I have given it only a cursory check, so I don't guarantee it works
correctly.
Depending on how your particular NetCDF interface works, you might need
to copy arrays from NetCDF to ensure they are genuine ndarray objects.
For plotting, you might want to try matplotlib. I think you will find
it easier to use than GMT, especially if you are accustomed to matlab.
Eric
Trond Kristiansen wrote:
Hi again.
I have attached the function that the FOR loop is part of as a python file.
What I am trying to do is to create a set of functions that will read the
output files (NetCDF) from running the ROMS model (ocean model). The output
file is organized in xi (x-direction), eta (y-direction), and s
(z-direction) where the s-values are vertical layers and not depth. This
particular function (z_slice) will find the closest upper and lower s-layer
for a given depth in meters (e.g. -100). Then values from the two selcted
layers will be interpolated to create a new layer at the selected depth
(-100). The problem is that the s-layers follow the bathymetry and a
particular s-layer will therefore sometimes be above and sometimes below the
selected depth that we want to interpolate to. That's why I need a quick
script that searches all of the layers and find the upper and lower layers
for a given depth value (which is negative). The z_r is a 3D array
(s,eta,xi) that is created using the netcdf module.
The main goal of these set of functions is to move away from using matlab,
but also to speed things up. The sliced data array will be plotted using GMT
or pyNGL.
Thanks for helping me. Cheers, Trond
'''
Suppose you have a 3-D array of (negative) depths
with indices level, x, and y. You also have a field,
say T, indexed the same way, and you want to interpolate
it to a given depth. Here is a quick example of how to
do it with no loops. This is not carefully tested, may
be buggy or incorrect, needs docstrings, etc.
It is intended to work for arrays of any dimension greater than
1 provided the first dimension is the one indexed by level.
'''
import numpy as np
def brackets(z, z0):
below = (z > z0).astype(np.int16).sum(axis=0)
above = below - 1
badmask = ((below == z.shape[0]) | (above < 0)).astype(np.bool)
above[badmask] = 0
below[badmask] = z.shape[0] - 1
return below, above, badmask
def interp(T, z, z0):
below, above, badmask = brackets(z, z0)
ind = np.indices(z.shape[1:])
ibelow = tuple([below] + list(ind))
iabove = tuple([above] + list(ind))
Tb = T[ibelow]
zb = z[ibelow]
dT = T[iabove] - Tb
dz = z[iabove] - zb
Ti = Tb + (dT/dz) * (z0 - zb)
return Ti, badmask
#Example use:
# Make depth levels varying linearly from the bottom to the surface:
nlevels = 5 # Use larger values to verify convergence.
z0 = np.linspace(0.0, -5000.0, num=nlevels)
z0.shape = nlevels,1,1
z1 = np.linspace(0.5, 1.0, num=4)
z1.shape = 1,4,1
z2 = np.linspace(0.1, 1.0, num=3)
z2.shape = 1,1,3
z = z0*z1*z2
# Temperature increases quadratically from 0 at -5000 m to 20 at the surface:
T = 20.0 * (1 + z/5000.0)**2
T1500, T1500badmask = interp(T, z, -1500)
# Now you might want to use a masked array to keep track of
# the valid values:
T1500m = np.ma.array(T1500, mask=T1500badmask)
print T1500m
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
