Griddata interpolation of data with gaps in Python

0
votes

I have issues interpolating my data into a grid using scipy interpolate griddata. My data is in 2D and has the form: X array, Y array and some intensity in Z. As these data are not spaced evenly, I want to reinterpolate them into a grid, in order to save images to work with later.

It works perpectly if I have one Z data block. My problem arises when I have discontinous Z data (like "islands" of data separated by zones were there is no data). In this case, interpolation also occurs in the gap between my data blocks.

Some illustration of the problem:

Image interpolation problem

Is there any way to interpolate my data properly?

Some part of the (very simple) code:

import numpy as np
from scipy import interpolate

X=data['X']
Y=data['Y']
Z=data['Z']

Xa,Ya=np.linspace(min_x, max_x, dim_x),np.linspace(min_y, max_y, dim_y) #dimensions of my grid, depends on the dataset I have
XX,YY=np.meshgrid(Xa,Ya) #creation of the grid
Zb = interpolate.griddata((X,Y), Z, (XX,YY), method='linear')

I have tried 'nearest' or 'cubic' but I does not work either...

EDIT: I added a text file with some example X, Y an Z I tried to interpolate, along with the new "Zb" interpolated data to show the problem. The same data was plotted in the image above. The data is available here: Gdrive

2
pythonnumpyscipydatasetinterpolation
Could you provide some sample data with an example of what you want to achieve? I guess you could either not compute the interpolation where you don't need it (when you generate the XX,YY points) or if you can afford the extra unuseful computation you could mask the empty regions back after the interpolation - filippo
I added a link to some .txt files with some data I want to plot along with my result after interpolation. I modified the example image to fit what I want to achieve with the data I provided. I don't know how to not compute the interpolation where I don't need it, as I don't have the coordinates of the empty regions... - Al Menia
Why not just set all the points outside of the spheres to black/zero after interpolation on the grid? Or, split the dataset in two, compute griddata on each side of the image with fill_value=0 and then recombine? - user228395
do you already have a mask that identifies interesting regions where you want to interpolate? if not maybe you should first think about how to obtain that, maybe some kind of clustering or even just morphology operations could be enough - filippo

2 Answers

0
votes

I had a similar problem. I solved this by computing differents grid (in your problem this should be one grid for each disk) and contour plot but with the same whole range of values and plotting them together. May be not 'elegant' but it is OK. Alain

0
votes

The solution I found consisted in abandonning the idea of interpolating my data. Instead, when reading my data, I defined a ROI (min_x,max_x,resolution_x and same for y) that I fill by iterating on my data, with various options (min, max, average...). It might not be the fastest solution, but it works and does not create artifacts!

import numpy as np

min_x, max_x, dim_x = (xxx, XXX,x) #limits on the dataset
min_y, max_y, dim_y = (yyy, YYY, y)

sorting_method="Max" #the method can be configured by the user

def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return idx

X=data['X'] #This dataset must be sorted so that only the data falling into the ROI is kept
Y=data['Y']
Z=data['Z']

pos_x=np.arange(min_x,max_x,res_x)
pos_y=np.arange(min_y,max_y,res_y)


ROI=np.zeros((dim_y,dim_x))

for values in X:
        loc_x=find_nearest(pos_x,values)
        loc_y=find_nearest(pos_y,Y[i])
        if sorting_method=='Max':
            ROI[loc_y,loc_x,0]=max(ROI[loc_y,loc_x,0],Z[i])

Zb=ROI[::-1,:]

Thank you everyone for their precious input!

Al

Example image on disks