Prev: newnarxsp multiple inputs
Next: Accessing Level 2 M File S Function Dialog Parameters
From: Chris on 7 Jul 2010 08:26
I am using a code that is used for optimal interpolation. My goal is to interpolate sediment diameters using this code. The test data the code comes with works fine however when I use my data it does not. Now the difference between the data that comes with is a average from the past being interpolated with new observations. I don't have averages from the past so what I am doing is I simply took the average of all my data (123 points) and used that as a background field to interpolate off of. Could some one take a quick look at the code and see if you can spot anything that would be producing the "NaN" in my output matrix called 'Inc'. Thanks

% Load coordinates and observation as Lon, Lat and observation
% Load coordinates and observation as lon, lat and To

% Create mesh of observation
meshfunction = TriScatteredInterp(Lon, Lat, observation);
[LON, LAT] = meshgrid(Lon, Lat);
Tb = meshfunction(LON, LAT);

Lat = Lat';
Lon = Lon';
lat = lat';
lon = lon';
To = To';
observation = observation';

%% Calculate forecast (Tf)
Tf=griddata(LON,LAT,Tb,lon,lat);
%% Define mapping parameters and
L=6; sig=0.1; nblocks=25;

%% Define blocks
dx=(Lon(2)-Lon(1))/2; xmin=min(Lon)-dx; xmax=max(Lon)+dx; dx=(xmax-xmin)/nblocks;
dy=(Lat(2)-Lat(1))/2; ymin=min(Lat)-dy; ymax=max(Lat)+dy; dy=(ymax-ymin)/nblocks;
[xal,yab]=meshgrid(xmin:dx:xmax-dx,ymin:dy:ymax-dy);
xal=xal(:); xar=xal+dx; xol=xal-L; xor=xar+L;
yab=yab(:); yat=yab+dy; yob=yab-L; yot=yat+L;

%% Calculate increment
Inc=nan*ones(size(LON));
for ib=1:nblocks^2
Io=find(lon>=xol(ib)&lon<=xor(ib)&lat>=yob(ib)&lat<=yot(ib));
Ia=find(LON>=xal(ib)&LON<=xar(ib)&LAT>=yab(ib)&LAT<=yat(ib));

Roo=gaspari_cohn(lon(Io),lat(Io),lon(Io),lat(Io),L);
Rao=gaspari_cohn(LON(Ia),LAT(Ia),lon(Io),lat(Io),L);

Inc(Ia)=Rao*inv(Roo+sig^2*eye(length(Io)))*(To(Io)-Tf(Io))';
end


------------- and the gaspari_cohn function is :

function rho=gaspari_cohn(x1,y1,x2,y2,L)

z1=x1(:) +1i*y1(:); n1=length(z1);
z2=x2(:)'+1i*y2(:)'; n2=length(z2);

rho=nan*ones(n1,n2);
for k=1:n2;
t=abs(z1(:)-z2(k));
phi=zeros(n1,1);
I=find(t/L<=0.5); T=t(I)/L; phi(I)=1-20/3*T.^2+5*T.^3+8*T.^4-8*T.^5;
I=find(t/L>=0.5&t/L<=1); T=t(I)/L; phi(I)=((8*T.^2+8*T-1).*(1-T).^4)/3./T;
rho(:,k)= phi./(1+(t/L).^2);
end
From: Steven Lord on 7 Jul 2010 09:49

"Chris " <chris.veinot(a)hotmail.com> wrote in message
news:i11rp6$ng4$1(a)fred.mathworks.com...
>I am using a code that is used for optimal interpolation. My goal is to
>interpolate sediment diameters using this code. The test data the code
>comes with works fine however when I use my data it does not. Now the
>difference between the data that comes with is a average from the past
>being interpolated with new observations. I don't have averages from the
>past so what I am doing is I simply took the average of all my data (123
>points) and used that as a background field to interpolate off of. Could
>some one take a quick look at the code and see if you can spot anything
>that would be producing the "NaN" in my output matrix called 'Inc'. Thanks

*snip*

> %% Calculate increment
> Inc=nan*ones(size(LON));

Well, this line looks suspicious :)

On a more serious note, the NaN function accepts a size input (and has for
quite a while) to generate an array of NaNs.

Inc = NaN(size(LON));

> for ib=1:nblocks^2
> Io=find(lon>=xol(ib)&lon<=xor(ib)&lat>=yob(ib)&lat<=yot(ib));

Looking below, you never actually use the _values_ of Io. You almost always
use Io as an index. In this case, you don't need to FIND; instead use Io to
perform logical indexing. The same holds for Ia.

> Ia=find(LON>=xal(ib)&LON<=xar(ib)&LAT>=yab(ib)&LAT<=yat(ib));
>
> Roo=gaspari_cohn(lon(Io),lat(Io),lon(Io),lat(Io),L);
> Rao=gaspari_cohn(LON(Ia),LAT(Ia),lon(Io),lat(Io),L);
>
> Inc(Ia)=Rao*inv(Roo+sig^2*eye(length(Io)))*(To(Io)-Tf(Io))';

First, DON'T INVERT. Use backslash instead, assuming Roo is a matrix of
size length(Io)-by-length(Io). [If you use the logical indexing suggestion
above, you'll need to change your EYE call to eye(nnz(Io)) or something
similar, to make an identity matrix whose number of rows and columns is the
number of true values in the logical vector Io.] This could, if
Roo+sig^2*eye(length(Io)) is ill-conditioned, return NaN values.

Second, what happens if there is a particular element of Inc that is never
"selected" as part of Ia? Wouldn't that element of Inc remain at its
preallocated value of NaN?

Set a breakpoint immediately after this loop and check Inc for NaN values,
then look at the corresponding elements of LON and LAT and your x*/y*
vectors to determine if the second case is the cause of the NaN values.

> end
>
>
> ------------- and the gaspari_cohn function is :
> function rho=gaspari_cohn(x1,y1,x2,y2,L)
>
> z1=x1(:) +1i*y1(:); n1=length(z1);
> z2=x2(:)'+1i*y2(:)'; n2=length(z2);
>
> rho=nan*ones(n1,n2);

See the preallocation tip above.

> for k=1:n2; t=abs(z1(:)-z2(k)); phi=zeros(n1,1);
> I=find(t/L<=0.5); T=t(I)/L;
> phi(I)=1-20/3*T.^2+5*T.^3+8*T.^4-8*T.^5;

See the logical indexing suggestion above.

*snip*

--
Steve Lord
slord(a)mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com


 | 
Pages: 1
Prev: newnarxsp multiple inputs
Next: Accessing Level 2 M File S Function Dialog Parameters