1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
|
**********************************************************************************************
*************************** Overview of the CDAT interface to natgrid ************************
**********************************************************************************************
---------------------------------------------------------------------------------------------
INTRODUCTION TO NGMATH
The ngmath library is a collection of interpolators and approximators for one-dimensional, two-dimensional
and three-dimensional data. The packages, which were obtained from NCAR, are:
natgrid -- a two-dimensional random data interpolation package based on Dave Watson's nngridr.
dsgrid -- a three-dimensional random data interpolator based on a simple inverse distance weighting
algorithm.
fitgrid -- an interpolation package for one-dimensional and two-dimensional gridded data based on
Alan Cline's Fitpack. Fitpack uses splines under tension to interpolate in one and two
dimensions.
csagrid -- an approximation package for one-dimensional, two-dimensional and three-dimensional random
data based on David Fulker's Splpack. csagrid uses cubic splines to calculate its
approximation function.
cssgrid -- an interpolation package for random data on the surface of a sphere based on the work of
Robert Renka. cssgrid uses cubic splines to calculate its interpolation function.
shgrid -- an interpolation package for random data in 3-space based on the work of Robert Renka.
shgrid uses a modified Shepard's algorithm to calculate its interpolation function.
COMPARISION OF NGMATH PACKAGES
Three-dimensional packages -- shgrid, csagrid and dsgrid.
shgrid is probably the package of choice for interpolation. It uses a least squares fit of biquadratics
to construct its interpolation function. The interpolation function will pass through the original data
points.
csagrid uses a least squares fit of cubic splines to calculate its approximation function: the calculated
surface will not necesarily pass through the original data points. The algorithm can become unstable in data
sparse regions.
dsgrid uses a weighted average algorithm and is stable in all cases, but the resultant interpolation is
not usually smooth and execution time is very slow. dsgrid is probably best used when csagrid and shgrid
fail or for comparative purposes.
Two-dimensional packages -- natgrid, fitgrid, csagrid and dsgrid.
natgrid is the package of choice in most cases. It implements a very stable algorithm and has parameters
for adjusting the smoothness of the output surface.
fitgrid offers user-settable parameters for specifiying derivatives along the boundary of the output grid
which are not available in natgrid.
csagrid produces an approximate two-dimensional surface which may be smoother than that produced by fitgrid
and natgrid.
dsgrid is not recommended for two-dimensional surfaces. natgrid is superior in all respects.
One-dimensional packages -- fitgrid and csagrid.
fitgrid is definitely the package of choice. It has many features not available in csagrid, such as
interpolating parametric curves, finding integrals, handling periodic functions, allowing smoothing that
varies from linear to a full cubic spline interpolation and specifying slopes at the end points.
Interpolation on a sphere -- cssgrid.
cssgrid is designed specifically for interpolating on a sphere. It uses cubic splines to calculate an
interpolation function.
NATGRID PACKAGE
natgrid implements a natural neighbor interpolation method. The input for the interpolation is a set
of randomly spaced two-dimensional coordinates with functional values at those coordinates; the output is a
set of interpolated values at coordinates in a user specified rectangular grid. The coordinates in the output
grid must be monotonic in each coordinate direction, but need not be evenly spaced. It is also possible to
interpolate at a single point.
natgrid uses a weighted average method that is much more sophisticated than the inverse distance weighted
average used by dsgrid. One distinguishing quality of natural neighbor interpolation is the way in which
a set of neighboring points (the natural neighbor) is selected to use for interpolating at a point. The
natural neighbor selection process avoids the problems common to methods based on choosing a fixed number
of neighboring points, or all points within a fixed distance. Another distinguishing quality of natural
neighbor interpolation is the way that the weights are calculated for the functional values at the natural
neighbor coordinates. These weights are based on proportionate area, rather than distances.
The method of finding the natural neighbors and calculating area-based weights to produce interpolated
values is called natural neighbor linear interpolation. This produces an interpolation surface that has a
continous slope at all points, except at the original input points. The result of natural neighbor linear
interpolation can be visualized as producing a snugly fit sheet stretched over all of the input points.
The interpolation method in natgrid also allows for natural neighbor linear interpolation augmented by
blending in gradient estimates. This is called natural neighbor nonlinear interpolation. It produces an
interpolation surface that has a continuous slope at all locations; two tautness parameters can be set by
the user to control the apparent smoothness of the output surface.
NATGRID CONTENTS
Access through Python to the natgrid package from NCAR's ngmath distribution is provided directly through the module
natgridmodule.so which was generated as a Python C language extension in order to export the natgrid functions
from the original C language library to Python.
REQUIRED FILE
natgridmodule.so -- the Python interface to the ngmath natgrid package.
USEFUL FILES
nat.py -- the object oriented interface including a general help package.
natgridtest.py -- the code to test nat.py and to write documentation.
USAGE
This module is designed to use in two ways. One is through the use of the object oriented interface to the underlying
functions. This approach is recommended for users not already familiar with the original natgrid distribtution because
it simplifies the calls to the routines. The other method uses the original functions calling them directly from Python.
------------------- OBJECT ORIENTED APPROACH ----------------
The nat module contains the Natgrid class and its single method, rgrd, which provides access to all the natgrid
functions. The object oriented approach has been organized as a two step process.
STEP 1.
To make an instance, r, type:
import nat
r = nat.Natgrid(xi, yi, xo, yo)
or
r = nat.Natgrid(xi, yi, xo, yo, listOutput = 'yes')
where xi, yi and xo, yo are the input and output grid coordinate arrays. The optional listOutput must
set to anything except 'no' if xo, yo are in list format as explained below. It is the responsibility
of the user to set listOutput if the output is in the list form.
The input grid must be organized in a list format always. The size of the xi array and the yi array are
necessarily equal. For example, if there are n randomly spaced input data points, there
are n values in xi and n values in yi.
There are two possible formats for the output grid. The output grid coordinate arrays may be a list like
the input array or it may be a rectangular grid. The choice between the two posibilities is made according
to requirements in subseqent calls to the method function. The first choice is required if the subsequent
call is to the single point mode interpolation. The list can have one or more points. Of course, the list
could describe a rectangular grid. For example, a rectangular grid with 10 x values and 20 y values can be
rewrtten in list form with 200 x value and 200 y values. However, this form requires calling the slower
single point interpolator. The second choice is most efficient for the basic interpolation to a rectangular
output grid. The output grid must be monotonic but need not be equally spced.
The grid coordinate arrays can be single precision (numpy.float32) or double precision (numpy.float64). The
decision on whether to call for a single or a double precision computation subsequently is made by looking at
the type of these arrays.
To look at the default settings for the control parameters and a brief description of thier properties, type
r.printDefaultParameterTable()
To change a setting type the new value. For example, to set igr to 1, type
r.igr = 1
To find a value without printing the table, type the name. For example, to exam the value of hor, type
r.hor
To check the settings type
r.printInstanceParameterTable() -- prints in tabular form the parameters used in subsequent calls to the method
function rgrd.
or
printStoredParameters() -- prints the parameters in memory which may differ from the above if the user
has made more than one instance of the Natgrid class.
STEP 2.
natgrid is restricted to two dimensions . Consequently, it is the user's responsibility to reduce the processing of
higher dimensional data to a sequence of calls using only two dimensional data.
The computations are divided into two groups depending on whether the output arrays are in list form or in rectilinear
grid form. If they are in list format the single point mode is called to interpolate to those individual points. This is
the only process possible. On the other hand, if the output goes to a rectangular grid there are more choices. In
addition to carrying out linear and nonlinear interpolations, it is possible to request aspects and slopes. The aspect
at a point on the interpolated surface is the direction of steepest descend. The slope is the value of the partial
derivative taken in the direction of the aspect. The slope is measured as an angle that is zero in a horizonal surface
and positive below the horizontal.
The following examples cover the basic computations. They start with a indication of the appropriate STEP 1.
Example 1: the basic natural neighbor linear interpolation
As STEP 1 make an instance, r, with:
import nat
r = nat.Natgrid(xi, yi, xo, yo)
where the xo, yo grid is rectilinear as explained above in STEP 1.
Then call the primary interpolation computation to regrid the input data, dataIn, on the grid (xi, yi) to
the output data, dataOut, on the grid (xo, yo), with
dataOut = r.rgrd( dataIn )
The computation is either single or double precision as determined by the precision submitted in the grid
description in STEP 1.
It is also possible to request a wrap in the input grid and the input data in the longitude direction, assumed
to be the yi grid coordinate, by adding a keyword as
dataOut = r.rgrd( dataIn, wrap = 'yes' )
Example 2: natural neighbor linear interpolation returning the aspect and the slope.
As STEP 1 make an instance, r, with:
import nat
r = nat.Natgrid(xi, yi, xo, yo)
where the xo, yo grid is rectilinear as explained above in STEP 1.
Then call the primary interpolation computation to regrid the input data, dataIn, on the grid (xi, yi) to
the output data, dataOut, on the grid (xo, yo), while asking for the aspect and the slope on this output grid, with
dataOut, a, s = r.rgrd( dataIn, aspectSlope = 'yes' )
where a is the aspect, the direction of the steepest descent in degrees measured from 'north' and s is the
slope in degrees measured from the horizontal. Necessarily, these are arrays aligned with the rectilinear
output grid, xo, yo.
The computation is either single or double precision as determined by the precision submitted in the grid
description in STEP 1.
It is also possible to request a wrap in the input grid and the input data in the longitude direction, assumed
to be the yi grid coordinate, by adding a keyword as
dataOut, a, s = r.rgrd( dataIn, aspectSlope = 'yes', wrap = 'yes' )
Example 3: the basic natural neighbor nonlinear interpolation
The procedure for the nonlinear interpolation differs from the linear case in the need to set the control
parameter igr. Follow Example 1 and insert the following statament after making the instance, r.
r.igr = 1
Example 4: natural neighbor nonlinear interpolation returning the aspect and the slope.
The procedure for the nonlinear interpolation differs from the linear case in the need to set the control
parameter igr. Follow Example 2 and insert the following statament after making the instance, r.
r.igr = 1
Example 5: single point mode natural neighbor linear interpolation
As STEP 1 make an instance, r, with:
import nat
r = nat.Natgrid(xi, yi, xo, yo, listOutput = 'yes')
where the xo, yo output grid is in the list form (not a rectangular output grid) as explained above in
STEP 1.
To call the single point mode interpolation computation to regrid the input data, dataIn, on the grid (xi, yi)
to the output data, dataOut, on the grid (xo, yo), type
dataOut = r.rgrd( dataIn )
The computation is either single or double precision as determined by the precision submitted in the grid
description in STEP 1. In the single point mode it is not possible to request the aspect and the slope.
Example 6: single point mode natural neighbor nonlinear interpolation
The procedure for the nonlinear interpolation differs from the linear case in the need to set the control
parameter igr. Follow Example 5 and insert the following statament after making the instance, r.
r.igr = 1
------------------- ORIGINAL FUNCTION APPROACH -----------------
The module natgridmodule.so exports the following functions to Python from the original ngmath C library:
Single precision procedures:
natgrids - primary function for gridding.
seti - set int parameter values.
geti - retrieve values for int parameters.
setr - set float parameter values.
getr - retrieve values for float parameters
setc - set char parameter values.
getc - retrieve values for char parameters.
getaspects - get aspect values, if calculated by setting sdi = 1.
getslopes - get slope values, if calculated by setting sdi = 1.
pntinits - initiate single point mode.
pnts - interpolate at a single point.
pntend - terminate single point mode.
Double precision procedures:
natgridd - primary function for gridding.
setrd - set float parameter values.
getrd - retrieve values for float parameters
getaspectd - get aspect values, if calculated by setting sdi = 1.
getsloped - get slope values, if calculated by setting sdi = 1.
pntinitd - initiate single point mode.
pntd - interpolate at a single point.
pntendd - terminate single point mode.
Information on the use of the routines is available by importing natgridmodule and printing the docstring
of interest. For example, documentation for the routine natgrids is obtained by typing
import natgridmodule
print natgridmodule.natgrids.__doc__
This same information is available in the help package.
A description of the control parameters is not in the natgridmodule documentation. It can be found by typing
import nat
nat.printParameterTable()
The documentation associated with the natgridmodule.so, such as the doctrings, describe the C code.
DOCUMENTATION
Documentation is provided through Python's docstrings, essentially Python style program
comments. A help package provides instructions on the use of the natgrid module. A table of contents
is printed to the screen by typing
nat.help()
after importing nat.
A hard copy of all the pertinent 'docstring' documentation written to the file natgridmodule.doc can
be produced by typing
nat.document()
As an alternate to using the help package, online documentation for the natgrids function, for example,
is available directly from the natgrids doctring by typing
import natgridmodule
print natgridmodule.natgrids.__doc__
TESTING
To run a test of the natgrid computations and to get a copy of this documentation, type
cdat natgridtest.py
--------------------------------------------------------------------------------------------------------------
HELP PACKAGE EXAMPLE
************************ Default Parameter Table **********************
-----------------------------------------------------------------------------------------------------
name type legal value default value description
---- ----- ------------ --------------- ------------------------------------------------------------
adf int 0 = no or 1 0 produce data file of algoritmic info for display? (see alg)
alg char any file nam "nnalg.dat" file name for algoritmic display tool (see adf)
asc int 0 = no or 1 1 is automatic scaling is allowed?
bI float >= 1. 1.5 tautness increasing effect of the gradients by increasing bI
bJ float >= 1. 7.0 tautness decreasing breadth of region affected by gradients
dup int 0 = yes or 1 1 are duplicate input coordinates are allowed?
ext int 0 = no or 1 1 is extrapolation allowed outside the convex hull?
hor float >= 0. -1.0 amount of horizontal overlap from outside current region
igr int 0 = no or 1 0 are gradients are to be computed?
magx float > 0. 1.0 scale factor for x coordinate values
magy float > 0. 1.0 scale factor for y coordinate values
magz float > 0. 1.0 scale factor for z coordinate values
non int 0 = yes or 1 0 are interpolated values are allowed to be negative?
nul float any float 0.0 value for points outside the convex hull if no extrapolation
rad int 0 = rad or 1 0 are slopes and aspects are returned in radians or degrees?
sdi int 0 = no or 1 0 are slopes and aspects to be computed?
upd int 0=N to S or 1 does output array from giving N to S or S to N?
ver float >= 0. -1.0 amount of vertical overlap from outside current region
xas float > 0. 0.0 scale used by automatic scaling of x in last interpolation
yas float > 0. 0.0 scale used by automatic scaling of y in last interpolation
zas float > 0. 0.0 scale used by automatic scaling of z in last interpolation
|
