Inverse distance weighting on scatter data sets in matlab

I want to use IDW interpolation technique on my data set.As usual consist on text files. Name Output_00.text to Output_23.text in a folder. Each text file consist on three columns. Latitude, Longitude and Temperature values. 3rd Column( Temperature column) contain -9999.000 value encoded as missing or NaN value. I want to interpolate these NaN/missing value in temperature through inverse distance weighting interpolation technique.

Here is what i am trying

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%INVERSE DISTANCE WEIGHT %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[Vint]=IDW(xc,yc,vc,x,y,e,r1,r2)
%%%%%%%%%%%%%%%%%
%%%INPUTS
%xc = stations x coordinates (columns) [vector]
%yc = stations y coordinates (rows) [vector]
%vc = variable values on the point [xc yc]
%x = interpolation points  x coordinates [vector]
%y = interpolation points y coordinates [vector]
%e = distance weight
%r1 --- 'fr' = fixed radius ;  'ng' = neighbours
%r2 --- radius lenght if r1 == 'fr' / number of neighbours if  r1 =='ng'
%%%OUTPUTS
%Vint --- Matrix [length(y),length(x)] with interpolated  variable values
%%%EXAMPLES
%%%--> V_spa=IDW(x1,y1,v1,x,y,-2,'ng',length(x1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Simone Fatichi -- [e-mail]
%   Copyright 2009
%   $Date: 2009/06/19 $ 
%   $Updated 2012/02/24 $ 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vint=zeros(length(y),length(x));
xc=reshape(xc,1,length(xc));
yc=reshape(yc,1,length(yc));
vc=reshape(vc,1,length(vc));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if  strcmp(r1,'fr')
    if  (r2<=0)
        disp('Error: Radius must be positive')
        return
    end
    for i=1:length(x)
        for j=1:length(y)
            D=[]; V=[]; wV =[]; vcc=[];
            D= sqrt((x(i)-xc).^2 +(y(j)-yc).^2);
            if min(D)==0
                disp('Error: One or more stations have the coordinates of an interpolation point')
                return
            end
            vcc=vc(D<r2); D=D(D<r2); 
            V = vcc.*(D.^e);
            wV = D.^e;
            if isempty(D)
                V=NaN;
            else
                V=sum(V)/sum(wV);
            end
            Vint(j,i)=V;
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
    if (r2 > length(vc)) || (r2<1)
        disp('Error: Number of neighbours not congruent with data')
        return
    end
    for i=1:length(x)
        for j=1:length(y)
            D=[]; V=[]; wV =[];vcc=[];
            D= sqrt((x(i)-xc).^2 +(y(j)-yc).^2);
            if min(D)==0
                disp('Error: One or more stations have the coordinates of an interpolation point')
                return
            end
            [D,I]=sort(D);
            vcc=vc(I);
            V = vcc(1:r2).*(D(1:r2).^e);
            wV = D(1:r2).^e;
            V=sum(V)/sum(wV);
            Vint(j,i)=V;
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return

I want to modify this code so that it read all text file step by step and interpolate NaN values with IDW interpolation technique and save each text file with IDW_00.text to IDW_23.text in the same folder.

I have attached data in text files. Please help and a lot of thanks for this kind assistance

[Muhammad Usman Saleem]

Please suggest me how can i set weight and search radius for interpolation NaN values in my data set. Problem is that NaN values appeared randomly and with out any sequence.I am also not know how many NaN values will be appear in text file. Guide me the best and convenient way of this solution.

[Muhammad Usman Saleem]

I am confusing in

%x = interpolation points  x coordinates [vector]
%y = interpolation points y coordinates [vector]
%e = distance weight
%r1 --- 'fr' = fixed radius ;  'ng' = neighbours
%r2 --- radius lenght if r1 == 'fr' / number of neighbours if  r1 =='ng'

How to set x and y coordinates of NaN values? what is r1 and r2?

Answer 1

[Ced]

I think you should ask yourself first why you want to use this method. What does the method do? What are it's advantages / disadvantages vs a typical cubic method? I would actually even start with a cubic spline, and compare. This should help to build your intuition.

To your questions:

a) How to set x and y coordinates of NaN values:

You already have the x and y coordinates of the NaN values, right? These are the latitude and longitude values. The problem is that you don't know the temperature. And that's what the interpolation is going to give you. So, in your algorithm:

xc, yc, vc: longitude, latitude and temperature of non-nan values.

x, y: longitude and latitude coordinates of your nan-values. This is where you want the algorithm to provide you with an interpolated guess.

b) What is r1 and r2?

From the code, we see the following: The IWD metric computes the distance between your desired interpolation points (in your case your nan values) and the known neighbours around them. The value of the neighbours will be weighted according to this distance. The distance used here is just the euclidean distance for e = 1. The type of norm (i.e. shape of your distance metric) is determined by e. I would start with e = 1.

r1: is used as a flag. There are two modes for this function. The meaning of r2 changes depending on the mode.

if r1 = 'fr', then r2 is the smallest weight that a point must have to be considered for the interpolation. It's the cutoff distance.

if r1 = 'ng', then r2 is the number of closest neighbours used for the interpolation.

Now for the actual values, you will have to tune them a bit to see what works best for you. Initially, you might start with r1 = 'fr', r2 = 10, e = 1.

[Muhammad Usman Saleem]

EDIT 2

I am trying this for one text file.

% Load the data
data = load('Output_00.txt');
% separate the data columns, just to make the code clear
Lat = data(:,1); % Column 1 is Latitude
Lon = data(:,2); % Column 2 is Longitude
Tmp = data(:,3); % Column 3 is Temperature
% Find the "good" data points
good_temp = find(Tmp > -9999.000);
% find the "bad" data points
bad_temp = find(Tmp == -9999.000);

now after that in the same folder i am calling this function as

Vint = IDW(Lat(good_temp),Lon(good_temp),Tmp(good_temp),Lat(bad_temp),Lon(bad_temp),1,'fr',10);

but error remain same

Error: File: IDW.m Line: 4 Column: 23
Unbalanced or unexpected parenthesis or bracket.

??

[Ced]

Just copy the whole "function [Vint] = IDW .... end" code (i.e. the actual function) into a new file and call it IDW.m

Then, from this "main" script, you call it as

Vint = IDW(Lat(good_temp),Lon(good_temp),Tmp(good_temp),Lat(bad_temp),Lon(bad_temp),1,'fr',10);

[Muhammad Usman Saleem]

After correcting the input arguments. I am getting that error?

>> Vint = IDW(Lat(good_temp),Lon(good_temp),Tmp(good_temp),Lat(bad_temp),Lon(bad_temp),1,'fr',10);
Error: One or more stations have the coordinates of an interpolation point

How to resolve this error?

[Muhammad Usman Saleem]

I think in my data set NaN values are in frequent and random in appearance. When i radius of 50 then also the error remains same. I think we should use portion of our concern from IDW function and interpolate NaN values in my datasets through IDW interpolation? But which portion will select?

[Ced]

The error means that some of the coordinates in xc, yc are identical to the ones in x, y.

[Muhammad Usman Saleem]

Thanks for your kind reply. Please check in my code above where i am producing mistake?

[Ced]

No time right now, sorry. I'll have a look at it tonight if you don't find the error in the meantime.

[Ced]

Hi

So.. the IDW function creates a whole grid of points, and none of them may match. Not sure what the idea is, but I am not familiar with the method, so ok. But since your "good" and "bad" grid are in the same domain, it fails. Also, you don't even want a full grid, you only want values for specific points.

I rewrote it into something else. Feel free to use it and see if it does something sensible ;). Be warned that this may be complete rubbish. Note that I only rewrote the first option, but the second one should be straightforward to adjust if you need it.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%INVERSE DISTANCE WEIGHT %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[Vint]=IDW(xc,yc,vc,x,y,e,r1,r2)
%%%%%%%%%%%%%%%%%
%%%INPUTS
%xc = stations x coordinates (columns) [vector]
%yc = stations y coordinates (rows) [vector]
%vc = variable values on the point [xc yc]
%x = interpolation points  x coordinates [vector]
%y = interpolation points y coordinates [vector]
%e = distance weight
%r1 --- 'fr' = fixed radius ;  
%r2 --- radius lenght if r1 == 'fr'
%%%OUTPUTS
%Vint --- Vector [length(1),1] with interpolated  variable values
%%%EXAMPLES
%%%--> V_spa=IDW(x1,y1,v1,x,y,1,'fr',10);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Adapted from Simone Fatichi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vint=zeros(length(x),1);
xc=xc(:);
yc=yc(:);
vc=vc(:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if  strcmp(r1,'fr')
    if  (r2<=0)
        disp('Error: Radius must be positive')
        return
    end
    for i=1:length(x)
            D = sqrt((x(i)-xc).^2 +(y(i)-yc).^2);
            if min(D)==0
                disp('Warning: One or more stations have the coordinates of an interpolation point')
                Vint(i) = vc(i);
                continue;
            end
            vcc=vc(D<r2); D=D(D<r2); 
            V = vcc.*(D.^e);
            wV = D.^e;
            if isempty(D)
                Vint(i)=NaN;
                disp('Warning: No point was close enough, setting Nan')
                continue
            else
                V=sum(V)/sum(wV);
            end
            Vint(i)=V;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

Let me know if it works.

Cheers

[Muhammad Usman Saleem]

Thanks for your assistance. This code is running fine. But its changing all the original data set. My all data sets has changes. Please read IDW where i takes this code. Also looks on the other ready made exiting IDW IDW IDW available on mathwork central

[Muhammad Usman Saleem]

From gIDW

I try this but both IDW and gIDW giving different results. Although running finely but NaN values appear very large as compared with neighbor. Please also i do not want to run a function rather a code simple?

function Vi = gIDW(Xc,Yc,Vc,Xi,Yi,w,r1,r2)
%%OUTPUTS
% Vi:   (mandatory) [PxQ]       gIDW interpolated values
%                                --> P=1, Q=1 yields interpolation at one
%                                    point
%                                --> P>1, Q=1 yields interpolation at a
%                                    vector of points
%                                --> P>1, Q>1 yields interpolation at a
%                                    (ir)regular grid of points
%%INPUTS
% Xc:   (mandatory) [Nx1]       x coordinates of known points
% Yc:   (mandatory) [Nx1]       y coordinates of known points
% Vc:   (mandatory) [Nx1]       known values at [Xc, Yc] locations
% Xi:   (mandatory) [PxQ]       x coordinates of points to be interpolated
% Yi:   (mandatory) [PxQ]       y coordinates of points to be interpolated
% w:    (mandatory) [scalar]    distance weight
%                                --> w<0, for Inverse Distance Weighted
%                                    interpolation [IDW]
%                                --> w=0, for Simple Moving Average (only
%                                    if neighorhood size is local and not
%                                    global) [SMA]
% r1:   (optional)  [string]    neighbourhood type
%                                --> 'n'     (default) number of neighbours
%                                --> 'r'     fixed radius length
% r2:   (optional)  [scalar]    neighbourhood size
%                                --> number of neighbours,  if r1=='n'
%                                    default is length(Xc)
%                                --> radius length,         if r1=='r'
%                                    default is largest distance between known points
%%SYNTAX
% ================== IDW ==================
% all inputs:
%   Vi = gIDW(Xc,Yc,Vc,Xi,Yi,-2,'n',30);
% 6 inputs:
%   Vi = gIDW(Xc,Yc,Vc,Xi,Yi,-2);
%       --> r1='n'; r2=length(Xc);
% 7 inputs:
%   Vi = gIDW(Xc,Yc,Vc,Xi,Yi,-2,'n');
%       --> r2=length(Xc);
%   Vi = gIDW(Xc,Yc,Vc,Xi,Yi,-2,'r');
%       --> r2=largest distance between know points [Xi,Yi] (see D1 calculation)
% ================== SMA ==================
%   Vi = gIDW(Xc,Yc,Vc,Xi,Yi,0,'n',10);
% ============== Spatial Map ==============
%   Vi = gIDW(Xc,Yc,Vc,Xi,Yi,-2,'n',10);
%       -with Xi and Yi 2D arrays of coordinates relative to an (ir)regular
%        grid.
%%EXAMPLES
% Interpolation at one point location:
%   Vi = gIDW([1:1:10]',[2:2:20]',rand(10,1)*100,5.5,11,-2,'n');
% =======================================
% Interpolation at a regular grid of unknown points:
%   XYc = [1:1:10]';
%   Vc = rand(10,1)*100;
%   Xi = rand(50,50)*10;
%   Yi = rand(50,50)*10;
%   [Xi,Yi] = meshgrid(XYc);
%   Vi = gIDW(XYc,XYc,Vc,Xi,Yi,-2,'r',3);
%   hold on
%   mapshow(Xi,Yi,Vi,'DisplayType','surface')
%   colormap gray
%   scatter(XYc,XYc,Vc,'filled','MarkerFaceColor','g','MarkerEdgeColor','y')
%   axis([0,11,0,11])
%   hold off
%%CODE
% check consistency of input parameters
if ~(length(Xc)-length(Yc)==0) || ~(length(Xc)-length(Vc)==0)
    error('varargin:chk',['Vectors Xc, Yc and Vc are incorrectly sized!'])
elseif ~(length(Xi)-length(Yi)==0)
    error('varargin:chk',['Vectors Xi and Yi are incorrectly sized!'])
elseif nargin < 6
    error('varargin:chk',['Uncorrect number of inputs!'])
end

% build input parameters
if nargin ~=8
    if nargin < 7   % default is 'n'
        r1 = 'n';
        r2 = length(Xc);
    elseif nargin==7 & r1=='n'
        r2 = length(Xc);
    elseif nargin==7 & r1=='r'  %for 'r' default is largest distance between know points
        [X1,X2] = meshgrid(Xc);
        [Y1,Y2] = meshgrid(Yc);
        D1 = sqrt((X1 - X2).^2 + (Y1 - Y2).^2);
        r2 = max(D1(:));     % largest distance between known points
        clear X1 X2 Y1 Y2 D1
    end
else
    switch r1 
        case {'r', 'n'}
            %nothing
        otherwise
            error('r1:chk',['Parameter r1 ("' r1 '") not properly defined!'])
    end
end

% initialize output
Vi = zeros(size(Xi,1),size(Xi,2));
D=[]; Vcc=[];

% fixed radius
if  strcmp(r1,'r')
    if  (r2<=0)
        error('r2:chk','Radius must be positive!')
        return
    end
    wb = waitbar(0,mfilename);
    for i=1:length(Xi(:))
        waitbar(i/length(Xi(:)))
%         if length(Xi(:))> 100, progress_bar(i, length(Xi(:)), mfilename); end
        D = sqrt((Xi(i)-Xc).^2 +(Yi(i)-Yc).^2);
        D = D(D<r2);
        Vcc = Vc(D<r2);
        if isempty(D)
            Vi(i) = NaN;
        else
            if sum(D==0)>0
                Vi(i) = Vcc(D==0);
            else
                Vi(i) = sum( Vcc.*(D.^w) ) / sum(D.^w);
            end
        end
    end
    close(wb)
% fixed neighbours number
elseif  strcmp(r1,'n')
    if (r2 > length(Vc)) || (r2<1)
        error('r2:chk','Number of neighbours not congruent with data')
        return
    end
    wb = waitbar(0,mfilename);
    for i=1:length(Xi(:))
        waitbar(i/length(Xi(:)))
%         if length(Xi(:))>100, progress_bar(i, length(Xi(:)), mfilename); end
        D = sqrt((Xi(i)-Xc).^2 +(Yi(i)-Yc).^2);
        [D,I] = sort(D);
        Vcc = Vc(I);
        if D(1) == 0
            Vi(i) = Vcc(1);
        else
            Vi(i) = sum( Vcc(1:r2).*(D(1:r2).^w) ) / sum(D(1:r2).^w);
        end
    end
    close(wb)
end

return

[Ced]

I'm pretty sure that

gIDW(xc,yx,vc,x,y,-2,'r',10)

and my version of

IDW(xc,yc,zc,x,y,-2,'fr',10)

should return the same answer. If not, take his version, mine was hacked together quickly. I guess that e = w > 0 does not make sense in this context. If you think that the returned values are not how you like them, then you need to tune the parameters.

Do you know that you can set breakpoints in matlab functions? Just click on the line number of the left in your editor. This will make a red dot appear (after you have saved). Next time you run the code, the code will stop there, allowing you to inspect the different variables etc.

That way, you can debug and inspect your code step by step.

Well, as I said in the beginning... you need to decide what you want to do. But if you want to use IDW, I would strongly suggest using functions. Pasting everything in your script is going to make things complete illegible.

[Muhammad Usman Saleem]

I have data of 24 such text files in 30 different folders. Running function and giving inputs each times will consume time. That why i am insisting over scripting. Tell me if short solution is available?

[Ced]

But that is exactly why you should use a function instead of a script. It's to keep things compact since you are performing the same operation on each file (and/or folder). What you could do is the following:

1. Write a script which has all the filenames and folders you want to load (e.g. in a cell of strings or something)

2. Write a loop which for each file, calls a function which extracts the "good" and "bad" samples, and then calls the IDW function to interpolate.

This way, all your files and folders are handled in a single file, and the whole loop will only have a couple of lines.

[Muhammad Usman Saleem]

no idea what you are saying?

[Ced]

I'm saying your code should look something like this:

clearvars

% 1. Have a cell which contains all the files you want to load

filenames = { 'file1.txt', 'file2.txt', 'file3.txt', ... };
N_files = length(filenames);

% 2. Now, for each file, load the data, and perform 
% whatever processing you need
data_final = cell(N_files,1); % cell to save the data

for i = 1:N_files
   % Load the data
   data = load(filenames{i});

   % separate the data columns, just to make the code clear
   Lat = data(:,1); % Column 1 is Latitude
   Lon = data(:,2); % Column 2 is Longitude
   Tmp = data(:,3); % Column 3 is Temperature
   % Find the "good" data points
   good_temp = find(Tmp > -9999.000);
   % find the "bad" data points
   bad_temp = find(Tmp == -9999.000);

   % perform the interpolation
   Tmp(bad_temp) = gIDW(Lat(good_temp), Lon(good_temp), ...
            Tmp(good_temp), Lat(bad_temp),Lon(bad_temp),-2,'r',10);
   % Save the data, rewrite to file, process... whatever you need
   data_final{i}.Lat = Lat;
   data_final{i}.Lon = Lon;
   data_final{i}.Tmp = Tmp;
end

As you can see, all your files and folders are easily handled in a loop. No need to copy the whole gIDW code in there.

[Muhammad Usman Saleem]

for exporting i am using

data1=[Lat';Lon';Tmp']';%here is 3*323 matrix of each file
 % Export.
     fId = fopen( sprintf( 'gIDW_interpolated%d.txt', i), 'w' ) ;
     fprintf( fId, '\t\t\r\n' ) ;
     fprintf( fId, '%.3f\t%.3f\t%.3f\r\n', data1(:) ) ;
     fclose( fId ) ;

Not sure whether data1(:) will be come in second last line or not??. I am getting 24 text files but latitude, longitude and temperature and randomly disturbed in each text file. See attachment

Getting that error

Undefined variable datal.

[Ced]

hm... the writing looks correct to me. I believe you have a typo somehwere. Note that the error says datal (with an L) and not data1 (with a 1). Note also that you still have all the "-9999.9999" in your file, which should not be there anymore! You are most likely writing the wrong matrix to file.

[Muhammad Usman Saleem]

I try this

clearvars

% 1. Have a cell which contains all the files you want to load
%pressures = [1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, 15, 10, 7, 5, 3, 2, 1.5, 1] ;
filenames = { 'Output_00.txt', 'Output_01.txt', 'Output_02.txt', 'Output_03.txt','Output_04.txt','Output_05.txt','Output_06.txt','Output_07.txt','Output_08.txt','Output_09.txt','Output_10.txt','Output_11.txt','Output_12.txt','Output_13.txt','Output_14.txt','Output_15.txt','Output_16.txt','Output_17.txt','Output_18.txt','Output_19.txt','Output_20.txt','Output_21.txt','Output_22.txt','Output_23.txt' };
N_files = length(filenames);

% 2. Now, for each file, load the data, and perform 
% whatever processing you need
data_final = cell(N_files,1); % cell to save the data

for i = 1:N_files
   % Load the data
   data = load(filenames{i});

     % separate the data columns, just to make the code clear
     Lat = data(:,1); % Column 1 is Latitude
     Lon = data(:,2); % Column 2 is Longitude
     Tmp = data(:,3); % Column 3 is Temperature
     % Find the "good" data points
     good_temp = find(Tmp > -9999.000);
     % find the "bad" data points
     bad_temp = find(Tmp == -9999.000);

     % perform the interpolation
     Tmp(bad_temp) = gIDW(Lat(good_temp), Lon(good_temp),Tmp(good_temp), Lat(bad_temp),Lon(bad_temp),-2,'r',10);
     % Save the data, rewrite to file, process... whatever you need
     data_final{i}.Lat = Lat;
     data_final{i}.Lon = Lon;
     data_final{i}.Tmp = Tmp;
  data1=[Lat';Lon';Tmp']';
     % Export.
       fId = fopen( sprintf( 'gIDW_inter%d.txt', i), 'w' ) ;
       fprintf( fId, '\t\t\r\n' ) ;
       fprintf( fId, '%.3f\t%.3f\t%.3f\r\n', data1) ;
       fclose( fId ) ;
   end

I got 23 text files but in each text file latitude is repeated in three columns. Then row of reparation of longitude come and repeat in three columns and in the last junk of interpolated temperature values appear in three columns. Tell me whether i am doing mistake? But when i call data1 in workspace it show me correct data in three columns lat,long, temp?

[Muhammad Usman Saleem]

i try also

fprintf( fId, '%.3f\t\n%.3f\t\n%.3f\r\n', data1) ;

but error remain same. I also change

data1=[Lat;Lon;Tmp];

but error remains same

[Muhammad Usman Saleem]

i just update the code with

data = [reshape(Lat.',[],1), reshape(Lon.',[],1), reshape(Tmp.',[],1)].' ;
   % Export.
   fId = fopen( sprintf( 'gIDW_inter%02d.txt', i-1), 'w' ) ;
     fprintf( fId, '\t\t\r\n' ) ;
     fprintf( fId, '%.3f\t%.3f\t%.3f\r\n', data(:)) ;
     fclose( fId ) ;

now this code storing in correct format. Please check line of data and tell me in explanations what this line is doing? Then i shall be very thankful to you

[Ced]

This is a very complicated way of transposing your matrix. You transpose it with .', then transpose is again with reshape, and then transpose it one last time with .'

What is the output if you replace the first line by

data = [ Lat Lon Tmp ]';

?

[Muhammad Usman Saleem]

thanks it is finely running.Please i have another query regarding to reading whole text file from a structure. Will you assist me in solution of this query? question is here

[Ced]

Nice! You can "accept" the answer to mark it as solved. I'll have a look at it tonight.