Download: https://osdyn.ifremer.fr/pyweb/notebooks/utils/interp/regrid_fortran_linear4dto1dxx.ipynb
fortran function linear4dto1dxx
[1]:
import numpy as N
from fortran._interp_ import linear4dto1dxx
note : Use of N.ogrid
[2]:
nti = 4
# inclusive
N.ogrid[0:nti-1:nti*1j]
[2]:
array([0., 1., 2., 3.])
[3]:
N.ogrid[0:nti+1:nti]
[3]:
array([0, 4])
[4]:
N.ogrid[0:nti-1:nti]
[4]:
array([0])
Basic initializations for the test
[5]:
# %% Bases
nex = 3
nexz = 2
nxi = 7
nyi = 6
nzi = 5
nti = 4
mv = 1.e20
no = 1
vfunc = lambda t, z, y, x: (1*x + 2.35*y + 3.24*z -0.65*t)
# Create t,z,y,x axes over a 4D grid
#1j way to use N.ogrid (inclusive)
tti, zzi, yyi, xxi = N.mgrid[0:nti-1:nti*1j, 0:nzi-1:nzi*1j,
0:nyi-1:nyi*1j, 0:nxi-1:nxi*1j]
print("Shapes of tti: {}, zzi: {}, yyi: {}, xxi: {}".format(tti.shape,zzi.shape,yyi.shape,xxi.shape))
# add one dimension
print("zzi.shape after N.mgrid {}".format(zzi.shape))
zzi = zzi[None]
print("zzi.shape after zzi[None] {}".format(zzi.shape))
zzi = N.repeat(zzi, nexz, axis=0)
print("zzi.shape after N.repeat {}".format(zzi.shape))
Shapes of tti: (4, 5, 6, 7), zzi: (4, 5, 6, 7), yyi: (4, 5, 6, 7), xxi: (4, 5, 6, 7)
zzi.shape after N.mgrid (4, 5, 6, 7)
zzi.shape after zzi[None] (1, 4, 5, 6, 7)
zzi.shape after N.repeat (2, 4, 5, 6, 7)
[6]:
nxi
[6]:
7
[7]:
xxi[0,0,0,:]
[7]:
array([0., 1., 2., 3., 4., 5., 6.])
[8]:
xxi[1,1,4,:]
[8]:
array([0., 1., 2., 3., 4., 5., 6.])
tti varies along 3th axis, same t values on z-y-x-dimensions
zzi varies along 2nd axis, same z values on t-y-x-dimensions
yyi varies along 1st axis, same y values on t-z-x-dimensions
xxi varies along 0th axis, same x values on t-z-y-dimensions
[9]:
N.testing.assert_equal(tti[3,...],tti[1,...]*3)
[10]:
N.testing.assert_equal(yyi[...,5,:],yyi[...,1,:]*5)
[11]:
N.testing.assert_array_equal(zzi[0,...],zzi[1,...]) # after repeat, i.e. extension to the left
Interpolation from pure 1D axes for t,z,y and x axes to a unique location
[12]:
# %% Pure 1D axes along each dimension
xi = xxi[0, 0, 0:1, :] # (nyix=1,nxi)
yi = yyi[0, 0, :, 0:1] # (nyi,nxiy=1)
zi = zzi[0:1, 0:1, :, 0:1, 0:1] # (nexz=1,ntiz=1,nzi,nyiz=1,nxiz=1)
ti = tti[:, 0, 0, 0] # (nti)
print("Size of xi: {}, yi: {}, zi: {}, ti: {},".format(xi.shape,yi.shape,zi.shape,ti.shape))
# 4D data field
vi = vfunc(tti, zzi, yyi, xxi)
print("Size of vi: {} and zzi: {}".format(vi.shape,zzi.shape))
# extend the first dimension of vi
vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,nti,nzi,nyi,nxi)
print("Size of vi after N.resize: {}".format(vi.shape))
# Get a unique random location
N.random.seed(0) #makes the random numbers predictable
xyztomin = -0.5
# get one number (no=1) between -0.5 and dim-1.5. saved into a numpy.array
xo = N.random.uniform(xyztomin, nxi-1.5, no)
yo = N.random.uniform(xyztomin, nyi-1.5, no)
zo = N.random.uniform(xyztomin, nzi-1.5, no)
to = N.random.uniform(xyztomin, nti-1.5, no)
# reference
vo_truth = N.ma.array(vfunc(to, zo, yo, xo))
# interpolation
mv = 1e20
vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)
print("vo_interp :",vo_interp)
#vo_interp = N.ma.masked_values(vo_interp), mv)
# None if correct
N.testing.assert_almost_equal(vo_interp[0], vo_truth)
Size of xi: (1, 7), yi: (6, 1), zi: (1, 1, 5, 1, 1), ti: (4,),
Size of vi: (2, 4, 5, 6, 7) and zzi: (2, 4, 5, 6, 7)
Size of vi after N.resize: (3, 4, 5, 6, 7)
vo_interp : [[15.47564723]
[15.47564723]
[15.47564723]]
[13]:
vo_interp
[13]:
array([[15.47564723],
[15.47564723],
[15.47564723]])
[14]:
vo_truth
[14]:
masked_array(data=[15.47564723],
mask=False,
fill_value=1e+20)
Important
Pourquoi 3 valeurs ??? Parce que 5ème dimension avant tzyx est égale à 3. Dimension non prise en compte dans interpolation
apparté
[15]:
vi.shape[1:]
[15]:
(4, 5, 6, 7)
[16]:
nex
[16]:
3
[17]:
(nex, 9)+vi.shape[1:]
[17]:
(3, 9, 4, 5, 6, 7)
Interpolation from a single point in space (x=cst, y=cst), z and t axes are pure 1D axes to a unique location
[18]:
# %% Single point in space
xi = xxi[0, 0, 0:1, :1] # (nyix=1,nxi)
yi = yyi[0, 0, :1, 0:1] # (nyi,nxiy=1)
zi = zzi[0:1, 0:1, :, 0:1, 0:1] # (nexz=1,ntiz=1,nzi,nyiz=1,nxiz=1)
ti = tti[:, 0, 0, 0] # (nti)
vi = vfunc(tti, zzi, yyi, xxi)[:, :, :, :1, :1]
vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,nti,nzi,1,1)
print("Size of xi: {}, yi: {}, zi: {}, ti: {},".format(xi.shape,yi.shape,zi.shape,ti.shape))
print("Size of vi: {} and zzi: {}".format(vi.shape,zzi.shape))
N.random.seed(0)
xyztomin = -0.5
xo = N.random.uniform(xyztomin, nxi-1.5, no)
yo = N.random.uniform(xyztomin, nyi-1.5, no)
zo = N.random.uniform(xyztomin, nzi-1.5, no)
to = N.random.uniform(xyztomin, nti-1.5, no)
vo_truth = N.ma.array(vfunc(to, zo, yi[0], xi[0]))
mv = 1e20
vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)
#vo_interp = N.ma.masked_values(vo_interp, mv)
print("vo_interp :",vo_interp)
N.testing.assert_almost_equal(vo_interp[0], vo_truth)
Size of xi: (1, 1), yi: (1, 1), zi: (1, 1, 5, 1, 1), ti: (4,),
Size of vi: (3, 4, 5, 1, 1) and zzi: (2, 4, 5, 6, 7)
vo_interp : [[5.45429115]
[5.45429115]
[5.45429115]]
[19]:
vo_truth
[19]:
masked_array(data=[5.45429115],
mask=False,
fill_value=1e+20)
[20]:
N.testing.assert_equal(vo_interp[0], vo_truth)
No interpolation in time but x, y and z axes are pure 1D axes
[21]:
# %% Constant time
xi = xxi[0, 0, 0:1, :] # (nyix=1,nxi)
yi = yyi[0, 0, :, 0:1] # (nyi,nxiy=1)
zi = zzi[0:1, 0:1, :, 0:1, 0:1] # (ntiz=1,nzi,nyiz=1,nxiz=1)
ti = tti[:1, 0, 0, 0] # (1)
vi = vfunc(tti, zzi, yyi, xxi)[:, :1, :, :, :]
vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,1,nzi,nyi,nxi)
N.random.seed(0)
xyztomin = -0.5
xo = N.random.uniform(xyztomin, nxi-1.5, no)
yo = N.random.uniform(xyztomin, nyi-1.5, no)
zo = N.random.uniform(xyztomin, nzi-1.5, no)
to = N.random.uniform(xyztomin, nti-1.5, no)
vo_truth = N.ma.array(vfunc(ti, zo, yo, xo))
mv = 1e20
vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)
#vo_interp = N.ma.masked_values(vo_interp, mv)
N.testing.assert_almost_equal(vo_interp[0], vo_truth)
[22]:
print("Size of xi: {}, yi: {}, zi: {}, ti: {}".format(xi.shape,yi.shape,zi.shape,ti.shape))
print("Size of vi: {} and zzi: {}".format(vi.shape,zzi.shape))
print("vo_truth : {}".format(vo_truth))
print("vo_interp :",vo_interp)
Size of xi: (1, 7), yi: (6, 1), zi: (1, 1, 5, 1, 1), ti: (1,)
Size of vi: (3, 1, 5, 6, 7) and zzi: (2, 4, 5, 6, 7)
vo_truth : [16.21316943]
vo_interp : [[16.21316943]
[16.21316943]
[16.21316943]]
[23]:
ti
[23]:
array([0.])
Interpolation from 1D x, y and t axes, z is a 5D array to a unique location
[24]:
# %% Variable depth with 1D X/Y + T
xi = xxi[0, 0, 0:1, :] # (nyix=1,nxi)
yi = yyi[0, 0, :, 0:1] # (nyi,nxiy=1)
zi = zzi[:, :, :, :, :] # (nexz,ntiz=nti,nzi,nyiz=nyi,nxiz=nxi)
ti = tti[:, 0, 0, 0] # (nti)
vi = vfunc(tti, zzi, yyi, xxi)
vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,nti,nzi,nyi,nxi)
N.random.seed(0)
xyztomin = -0.5
xo = N.random.uniform(xyztomin, nxi-1.5, no)
yo = N.random.uniform(xyztomin, nyi-1.5, no)
zo = N.random.uniform(xyztomin, nzi-1.5, no)
to = N.random.uniform(xyztomin, nti-1.5, no)
vo_truth = N.ma.array(vfunc(to, zo, yo, xo))
mv = 1e20
vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)
#vo_interp = N.ma.masked_values(vo_interp, mv)
N.testing.assert_almost_equal(vo_interp[0], vo_truth)
[25]:
print("Size of xi: {}, yi: {}, zi: {}, ti: {}".format(xi.shape,yi.shape,zi.shape,ti.shape))
print("Size of vi: {} and zzi: {}".format(vi.shape,zzi.shape))
print("vo_truth : {}".format(vo_truth))
print("vo_interp :",vo_interp)
Size of xi: (1, 7), yi: (6, 1), zi: (2, 4, 5, 6, 7), ti: (4,)
Size of vi: (3, 4, 5, 6, 7) and zzi: (2, 4, 5, 6, 7)
vo_truth : [15.47564723]
vo_interp : [[15.47564723]
[15.47564723]
[ nan]]
[26]:
xo,yo,zo,to
[26]:
(array([2.79288102]),
array([3.07594683]),
array([1.9110535]),
array([1.13464955]))
[27]:
ti
[27]:
array([0., 1., 2., 3.])
[28]:
print("Size of xi: {}, yi: {}, zi: {}, ti: {}".format(xi.shape,yi.shape,zi.shape,ti.shape))
print("Size of vi: {} and zzi: {}".format(vi.shape,zzi.shape))
print("vo_truth : {}".format(vo_truth))
print("vo_interp :",vo_interp)
Size of xi: (1, 7), yi: (6, 1), zi: (2, 4, 5, 6, 7), ti: (4,)
Size of vi: (3, 4, 5, 6, 7) and zzi: (2, 4, 5, 6, 7)
vo_truth : [15.47564723]
vo_interp : [[15.47564723]
[15.47564723]
[ nan]]
Error
Pourquoi nan pour dernier élément ???
Interpolation from 2D X/Y - no other axes (pure curvilinear) - to a unique location
[29]:
# %% 2D X/Y with no other axes (pure curvilinear)
xi = xxi[0, 0, ...] # (nyix=nyi,nxi)
yi = yyi[0, 0, ...]# (nyi,nxiy=nxi)
zi = zzi[0:1, 0:1, 0:1, 0:1, 0:1] # (nexz=1,ntiz=1,1,nyiz=1,nxiz=1)
ti = tti[:1, 0, 0, 0] # (1)
vi = vfunc(tti, zzi, yyi, xxi)[:, :1, :1, :, :]
vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,1,1,nyi,nxi)
N.random.seed(0)
xyztomin = -0.5
xo = N.random.uniform(xyztomin, nxi-1.5, no)
yo = N.random.uniform(xyztomin, nyi-1.5, no)
zo = N.random.uniform(xyztomin, nzi-1.5, no)
to = N.random.uniform(xyztomin, nti-1.5, no)
vo_truth = N.ma.array(vfunc(ti, zi.ravel()[0], yo, xo))
mv = 1e20
# xi and yi are 2D
vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)
# xi and yi are vectors - rectilinear grid
vo_interp_rect = linear4dto1dxx(xi[:1],yi[:, :1],zi,ti,vi,xo,yo,zo,to)
N.testing.assert_almost_equal(vo_interp[0], vo_truth)
N.testing.assert_almost_equal(vo_interp_rect[0], vo_truth)
[30]:
print("Size of xi: {}, yi: {}, zi: {}, ti: {}".format(xi.shape,yi.shape,zi.shape,ti.shape))
print("Size of vi: {} and zzi: {}".format(vi.shape,zzi.shape))
print("vo_truth : {}".format(vo_truth))
print("vo_interp :",vo_interp)
Size of xi: (6, 7), yi: (6, 7), zi: (1, 1, 1, 1, 1), ti: (1,)
Size of vi: (3, 1, 1, 6, 7) and zzi: (2, 4, 5, 6, 7)
vo_truth : [10.02135608]
vo_interp : [[10.02135608]
[10.02135608]
[10.02135608]]
Interpolation from and to the same coordinates
[31]:
# %% Same coordinates
xi = xxi[0, 0, 0:1, :] # (nyix=1,nxi)
yi = yyi[0, 0, :, 0:1] # (nyi,nxiy=1)
zi = zzi[0:1, 0:1, :, 0:1, 0:1] # (nexz=1,ntiz=1,nzi,nyiz=1,nxiz=1)
ti = tti[:, 0, 0, 0] # (nti)
vi = vfunc(tti, zzi, yyi, xxi)
vi = N.resize(vi, (nex, )+vi.shape[1:]) # (nex,nti,nzi,nyi,nxi)
xyzo = N.meshgrid(xi, yi, zi, ti, indexing='ij')
xo = xyzo[0].ravel()
yo = xyzo[1].ravel()
zo = xyzo[2].ravel()
to = xyzo[3].ravel()
vo_truth = N.ma.array(vfunc(to, zo, yo, xo))
mv = 1e20
vo_interp = linear4dto1dxx(xi,yi,zi,ti,vi,xo,yo,zo,to)
#vo_interp = N.ma.masked_values(vo_interp, mv)
N.testing.assert_almost_equal(vo_interp[0], vo_truth)
[32]:
print("Size of xi: {}, yi: {}, zi: {}, ti: {}".format(xi.shape,yi.shape,zi.shape,ti.shape))
print("Size of vi: {} and zzi: {}".format(vi.shape,zzi.shape))
print("vo_truth : {}".format(vo_truth.shape))
print("vo_interp :",vo_interp.shape)
print('6*7*5*4',6*7*5*4)
Size of xi: (1, 7), yi: (6, 1), zi: (1, 1, 5, 1, 1), ti: (4,)
Size of vi: (3, 4, 5, 6, 7) and zzi: (2, 4, 5, 6, 7)
vo_truth : (840,)
vo_interp : (3, 840)
6*7*5*4 840