Inverse distance weighting on scatter data sets in matlab

Question

Muhammad Usman Saleem 2016-3-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/273298-inverse-distance-weighting-on-scatter-data-sets-in-matlab

评论： Ced 2016-3-21

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 -- simonef@dicea.unifi.it
%   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

2 个评论
显示无隐藏无

Muhammad Usman Saleem 2016-3-13

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 2016-3-15

在 MATLAB Online 中打开

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 2016-3-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/273298-inverse-distance-weighting-on-scatter-data-sets-in-matlab#answer_213813

编辑：Ced 2016-3-16

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.

23 个评论
显示 21更早的评论隐藏 21更早的评论

Ced 2016-3-17

编辑：Ced 2016-3-17

在 MATLAB Online 中打开

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 2016-3-17

在 MATLAB Online 中打开

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 2016-3-19

在 MATLAB Online 中打开

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 2016-3-20

编辑：Muhammad Usman Saleem 2016-3-20

在 MATLAB Online 中打开

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 2016-3-21

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 2016-3-21