2
$\begingroup$

I have three-dimensional radar reflectivity data obtained as voxels (scans, rays, altitudes). The data has been sampled at irregular spacings and I want to convert this into a regular grid. In addition, I have 3-D arrays for the latitude, longitude, and altitude with the same shape(scans, rays, altitudes). Python is my programming language of choice. Is using scipy's RegularGridInterpolator the best way to obtain a regular grid?

A colleague suggested that I first start with a 2-D interpolation (scans and rays for every altitude) and then proceed to do a 3-D interpolation. A bulk of the reflectivity data contains zeros(zero reflectivity) and NaNs(no data available) and whether any filter would be applicable to remove the zeros/NaNs would be nice to know. Any recommendations and references on this topic will be appreciated.

So for a sample data set - I have total of 22163680 points and out of which 266111 are finite values(greater than zero). The rest are zeros and NaNs.

$\endgroup$
  • 2
    $\begingroup$ So, to clarify, you have data obtained in 3D at scattered point locations and wonder if there is an advantage of doing 2-D interpolation (say in XY planes) and then interpolating this data on a regular grid 2-D slices to full 3D? Together with probably having some filters in-between the interpolations? $\endgroup$ – Anton Menshov Apr 16 '18 at 16:41
  • $\begingroup$ @AntonMenshov - I couldn't have said it better. Yes that is what is exactly what I am asking $\endgroup$ – gansub Apr 16 '18 at 16:42
  • 1
    $\begingroup$ Interpolation from unstructured data sets is not a new topic. What have you already done, read, and tried? $\endgroup$ – Wolfgang Bangerth Apr 16 '18 at 16:46
  • $\begingroup$ @WolfgangBangerth - I have tried the basic scipy examples for RegularGridInterpolator found in the internet and regular 1D interpolation. But these were not sparse data sets. Here the bulk of the data is zero or NA $\endgroup$ – gansub Apr 16 '18 at 16:48
  • 1
    $\begingroup$ But there is a difference between data that is not available and data that is zero. If you have data on a regular $(i,j,k)$ grid where a lot of data points are zero is still dense, even if you don't actually store the zeros. The point is that you can still query at these points and get a valid answer -- the data is not sparse. Only if the data is genuinely not available is the data set sparse -- and in that case, you'll have to use insight into the data to say what is supposed to happen between data points. $\endgroup$ – Wolfgang Bangerth Apr 17 '18 at 0:03
1
$\begingroup$

For scattered data at points with no structure, try inverse-distance-weighted-idw-interpolation-with-python. This combines scipy's fast K-d_trees with inverse-distance aka radial kernels:

$\qquad interpol( x ) = \sum w_i \, f_i$, sum over say 10 points nearest $x$
$\qquad \qquad \qquad \qquad w_i = {1 \over |x - x_i|^p}$ normalized so $\sum w_i = 1$ .
It's only 40 lines of code, so easy to understand and modify. Runtimes on a stock iMac:
0.2 sec to build a Kd tree for 170k points in 3d
25 sec to query 100k points.

For rectangular grids aka box grids in 2d, 3d and up, data at e.g. $m \times n \times k$ points $(x_i, y_j, z_k)$ with no NAs, see Intergrid, a wrapper for scipy.ndimage.map_coordinates.

(Needless to say, it's a good idea to check any interpolation. An easy method is to split the data into say 90 % $train$ + 10 % $test$ points, and compare $interpolate( \, train \to test \, )$ with $true( \, test \, )$ .)

See also:
gis.stackexchange.com/questions/tagged/interpolation+python on gis.stack
how-to-interpolate-scattered-data-to-a-regular-grid-in-python on earthscience.stack
pyresample, pyproj .

$\endgroup$
  • $\begingroup$ No. Zeros mean value of zero reflectivity. NA is different. It means data was not available for some reason. $\endgroup$ – gansub May 6 '18 at 9:00
  • $\begingroup$ @gansub, both 0 s and NAs drop out of the sum -- you build a K-d tree and interpolate from the 266111 data points. (Please comment on how this worked.) $\endgroup$ – denis May 7 '18 at 8:11

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.