Subversion Repositories lagranto.arpege

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
""" interpolate data given on an Nd rectangular grid, uniform or non-uniform.
2
 
3
documentation and original website:
4
https://denis-bz.github.io/docs/intergrid.html
5
 
6
Purpose: extend the fast N-dimensional interpolator
7
`scipy.ndimage.map_coordinates` to non-uniform grids, using `np.interp`.
8
 
9
Background: please look at
10
http://en.wikipedia.org/wiki/Bilinear_interpolation
11
http://stackoverflow.com/questions/6238250/multivariate-spline-interpolation-in-python-scipy
12
http://docs.scipy.org/doc/scipy-dev/reference/generated/scipy.ndimage.interpolation.map_coordinates.html
13
 
14
Example
15
-------
16
Say we have rainfall on a 4 x 5 grid of rectangles, lat 52 .. 55 x lon -10 .. -6,
17
and want to interpolate (estimate) rainfall at 1000 query points
18
in between the grid points.
19
 
20
        # define the grid --
21
    griddata = np.loadtxt(...)  # griddata.shape == (4, 5)
22
    lo = np.array([ 52, -10 ])  # lowest lat, lowest lon
23
    hi = np.array([ 55, -6 ])   # highest lat, highest lon
24
 
25
        # set up an interpolator function "interfunc()" with class Intergrid --
26
    interfunc = Intergrid( griddata, lo=lo, hi=hi )
27
 
28
        # generate 1000 random query points, lo <= [lat, lon] <= hi --
29
    query_points = lo + np.random.uniform( size=(1000, 2) ) * (hi - lo)
30
 
31
        # get rainfall at the 1000 query points --
32
    query_values = interfunc( query_points )  # -> 1000 values
33
 
34
What this does:
35
    for each [lat, lon] in query_points:
36
        1) find the square of griddata it's in,
37
            e.g. [52.5, -8.1] -> [0, 3] [0, 4] [1, 4] [1, 3]
38
        2) do bilinear (multilinear) interpolation in that square,
39
            using `scipy.ndimage.map_coordinates` .
40
Check:
41
    interfunc( lo ) -> griddata[0, 0],
42
    interfunc( hi ) -> griddata[-1, -1] i.e. griddata[3, 4]
43
 
44
Parameters
45
----------
46
    griddata: numpy array_like, 2d 3d 4d ...
47
    lo, hi: user coordinates of the corners of griddata, 1d array-like, lo < hi
48
    maps: a list of `dim` descriptors of piecewise-linear or nonlinear maps,
49
        e.g. [[50, 52, 62, 63], None]  # uniformize lat, linear lon
50
    copy: make a copy of query_points, default True;
51
        copy=False overwrites query_points, runs in less memory
52
    verbose: default 1: print a 1-line summary for each call, with run time
53
    order=1: see `map_coordinates`
54
    prefilter: 0 or False, the default: smoothing B-spline
55
              1 or True: exact-fit interpolating spline (IIR, not C-R)
56
              1/3: Mitchell-Netravali spline, 1/3 B + 2/3 fit
57
        (prefilter is only for order > 1, since order = 1 interpolates)
58
 
59
Non-uniform rectangular grids
60
-----------------------------
61
What if our griddata above is at non-uniformly-spaced latitudes,
62
say [50, 52, 62, 63] ?  `Intergrid` can "uniformize" these
63
before interpolation, like this:
64
 
65
    lo = np.array([ 50, -10 ])
66
    hi = np.array([ 60, -6 ])
67
    maps = [[50, 52, 62, 63], None]  # uniformize lat, linear lon
68
    interfunc = Intergrid( griddata, lo=lo, hi=hi, maps=maps )
69
 
70
This will map (transform, stretch, warp) the lats in query_points column 0
71
to array coordinates in the range 0 .. 3, using `np.interp` to do
72
piecewise-linear (PWL) mapping:
73
    50  51  52  53  54  55  56  57  58  59  60  61  62  63  # lo[0] .. hi[0]
74
 
75
 
76
`maps[1] None` says to map the lons in query_points column 1 linearly:
77
    -10  -9  -8  -7  -6  # lo[1] .. hi[1]
78
 
79
 
80
More doc: https://denis-bz.github.com/docs/intergrid.html
81
 
82
"""
83
# split class Gridmap ?
84
 
85
from time import time
86
# warnings
87
import numpy as np
88
from scipy.ndimage import map_coordinates, spline_filter
89
import collections
90
 
91
__version__ = "2013-07-01 jul denis"
92
__author_email__ = "denis-bz-py@t-online.de"  # comments welcome, testcases most welcome
93
 
94
#...............................................................................
95
class Intergrid:
96
    __doc__ = globals()["__doc__"]
97
 
98
    def __init__( self, griddata, lo, hi, maps=[], copy=True, verbose=1,
99
            order=1, prefilter=False ):
100
        griddata = np.asanyarray( griddata )
101
        dim = griddata.ndim  # - (griddata.shape[-1] == 1)  # ??
102
        assert dim >= 2, griddata.shape
103
        self.dim = dim
104
        if np.isscalar(lo):
105
            lo *= np.ones(dim)
106
        if np.isscalar(hi):
107
            hi *= np.ones(dim)
108
        assert lo.shape == (dim,), lo.shape
109
        assert hi.shape == (dim,), hi.shape
110
        self.loclip = lo = np.asarray_chkfinite( lo ).copy()
111
        self.hiclip = hi = np.asarray_chkfinite( hi ).copy()
112
        self.copy = copy
113
        self.verbose = verbose
114
        self.order = order
115
        if order > 1  and 0 < prefilter < 1:  # 1/3: Mitchell-Netravali = 1/3 B + 2/3 fit
116
            exactfit = spline_filter( griddata )  # see Unser
117
            griddata += prefilter * (exactfit - griddata)
118
            prefilter = False
119
        self.griddata = griddata
120
        self.prefilter = (prefilter == True)
121
 
122
        self.nmap = 0
123
        if maps != []:  # sanity check maps --
124
            assert len(maps) == dim, "maps must have len %d, not %d" % (
125
                dim, len(maps))
126
            for j, (map, n, l, h) in enumerate( zip( maps, griddata.shape, lo, hi )):
127
                if map is None  or  isintance(map, collections.Callable):
128
                    lo[j], hi[j] = 0, 1
129
                    continue
130
                maps[j] = map = np.asanyarray(map)
131
                self.nmap += 1
132
                assert len(map) == n, "maps[%d] must have len %d, not %d" % (
133
                    j, n, len(map) )
134
                mlo, mhi = map.min(), map.max()
135
                if not (l <= mlo <= mhi <= h):
136
                    print("Warning: Intergrid maps[%d] min %.3g max %.3g " \
137
                        "are outside lo %.3g hi %.3g" % (
138
                        j, mlo, mhi, l, h ))
139
        self.maps = maps
140
        self._lo = lo  # caller's lo for linear, 0 nonlinear maps
141
        shape_1 = np.array( self.griddata.shape, float ) - 1
142
        self._linearmap = shape_1 / np.where( hi > lo, hi - lo, np.inf )  # 25jun
143
 
144
#...............................................................................
145
    def __call__( self, X, out=None ):
146
        """ query_values = Intergrid(...) ( query_points npt x dim )
147
        """
148
        X = np.asanyarray(X)
149
        assert X.shape[-1] == self.dim, ("the query array must have %d columns, "
150
                "but its shape is %s" % (self.dim, X.shape) )
151
        Xdim = X.ndim
152
        if Xdim == 1:
153
            X = np.asarray([X])  # in a single point -> out scalar
154
        if self.copy:
155
            X = X.copy()
156
        assert X.ndim == 2, X.shape
157
        npt = X.shape[0]
158
        if out is None:
159
            out = np.empty( npt, dtype=self.griddata.dtype )
160
        t0 = time()
161
        np.clip( X, self.loclip, self.hiclip, out=X )  # inplace
162
            # X nonlinear maps inplace --
163
        for j, map in enumerate(self.maps):
164
            if map is None:
165
                continue
166
            if isinstance(map, collections.Callable):
167
                X[:,j] = map( X[:,j] )  # clip again ?
168
            else:
169
                # PWL e.g. [50 52 62 63] -> [0 1 2 3] --
170
                X[:,j] = np.interp( X[:,j], map, np.arange(len(map)) )
171
            # linear map the rest, inplace (affine_transform ?)
172
        if self.nmap < self.dim:
173
            X -= self._lo
174
            X *= self._linearmap  # (griddata.shape - 1) / (hi - lo)
175
        ## print "test X:", X[0]
176
 
177
#...............................................................................
178
        map_coordinates( self.griddata, X.T,
179
            order=self.order, prefilter=self.prefilter,
180
            mode="nearest",  # outside -> edge
181
                # test: mode="constant", cval=np.NaN,
182
            output=out )
183
        if self.verbose:
184
            print("Intergrid: %.3g msec  %d points in a %s grid  %d maps  order %d" % (
185
                (time() - t0) * 1000, npt, self.griddata.shape, self.nmap, self.order ))
186
        return out if Xdim == 2  else out[0]
187
 
188
    at = __call__
189
 
190
# end intergrid.py