Inference of Model

Load and import dataset

clear all; load('data/Density_inference.mat','data'); load('coastlines.mat'); addpath('functions/')
Gneiting = @(dist,time,range_dist,range_time,delta,gamma,beta) 1./( (time./range_time).^(2.*delta) +1 ) .* exp(-(dist./range_dist).^(2.*gamma)./((time./range_time).^(2.*delta) +1).^(beta.*gamma) );

Histogram relatively close to normal distribution -> no transofmration used

figure('position',[0 0 1000 400]);

subplot(1,2,1); histogram(data.us);

subplot(1,2,2); histogram(data.vs);

No clear non-stionarinarity with NNT -> no trend used

figure('position',[0 0 1000 400]);

subplot(2,1,1);plot(data.NNT,data.us,'.k'); ylabel('E-W speed [m/s]');

subplot(2,1,2);plot(data.NNT,data.vs,'.k');ylabel('N-S speed [m/s]');

xlabel('Normalized Night Time (NNT) (-1: sunset | 1: sunrise)');

Average Daily

fmulti.t = data.day;

fmulti.r = 1:data.nrad;

[fmulti.R,fmulti.T]=meshgrid(fmulti.r,fmulti.t);

G = findgroups(data.day_id);

fmulti.Mus = splitapply(@nanmean,data.us,G);

fmulti.Mvs = splitapply(@nanmean,data.vs,G);

fmulti.isnan = isnan(fmulti.Mus);

figure('position',[0 0 1000 400]); hold on;

subplot(2,1,1); imagesc(fmulti.Mus')

subplot(2,1,2); imagesc(fmulti.Mvs')

Block décomposition

fmulti.b = nan(size(fmulti.t));
for i_b=1:numel(data.block.date)
    if i_b==numel(data.block.date)
        fmulti.b(fmulti.t>=data.block.date(i_b) | fmulti.t<data.block.date(1)) = i_b;
    else
        fmulti.b(fmulti.t>=data.block.date(i_b) & fmulti.t<data.block.date(i_b+1)) = i_b;
    end
end
fmulti.B = repmat(fmulti.b,1,data.nrad);
data.day_b = fmulti.b;

fmulti_poly_t_degree=0;
fmulti_poly_d_degree=1;
X = trend(fmulti.T,data.lat(fmulti.R),data.lon(fmulti.R),fmulti_poly_t_degree,fmulti_poly_d_degree);
fmulti.f_trend = @(t,lat,lon,beta) reshape(trend(t(:),lat(:),lon(:),fmulti_poly_t_degree,fmulti_poly_d_degree) * beta(:), size(t));
fmulti.Mus_n = nan(size(fmulti.Mus)); fmulti.Mvs_n = fmulti.Mus_n;
fmulti.beta_u = nan(fmulti_poly_t_degree + fmulti_poly_d_degree*2 +1 ,numel(data.block.date));
fmulti.beta_v = fmulti.beta_u;

We loop over all block...

for i_b=1:numel(data.block.date)
    id = find(~fmulti.isnan & fmulti.B==i_b);
    fmulti.beta_u(:,i_b) = ( X(id,:)' * X(id,:) ) \ (X(id,:)' * fmulti.Mus(id) );
    fmulti.beta_v(:,i_b) = ( X(id,:)' * X(id,:) ) \ (X(id,:)' * fmulti.Mvs(id) );
    fmulti.Mus_n(id) = fmulti.Mus(id)-fmulti.f_trend(fmulti.T(id),data.lat(fmulti.R(id)),data.lon(fmulti.R(id)),fmulti.beta_u(:,i_b));
    fmulti.Mvs_n(id) = fmulti.Mvs(id)-fmulti.f_trend(fmulti.T(id),data.lat(fmulti.R(id)),data.lon(fmulti.R(id)),fmulti.beta_v(:,i_b));
end

figure('position',[0 0 1000 400]);

subplot(2,1,1); hold on

for i_b=1:numel(data.block.date)-1

fill([data.block.date(i_b) data.block.date(i_b) data.block.date(i_b+1) data.block.date(i_b+1)],[25 -25 -25 25],data.block.col(i_b,:))

end

fill([min(data.day) min(data.day) data.block.date(1) data.block.date(1)],[25 -25 -25 25],data.block.col(4,:))

fill([data.block.date(end) data.block.date(end) max(data.day) max(data.day) ],[25 -25 -25 25],data.block.col(4,:))

plot(fmulti.T(:),fmulti.Mus(:),'.k');

for i_b=1:numel(data.block.date)

id = find(~fmulti.isnan & fmulti.B==i_b);

t_tmp = sort(unique(fmulti.T(id)));

plot(t_tmp, fmulti.f_trend(t_tmp,mean(data.lat(fmulti.R(id)))*ones(numel(t_tmp),1),mean(data.lon(fmulti.R(id)))*ones(numel(t_tmp),1), fmulti.beta_u(:,i_b)),'-r','Linewidth',2)

end

ylabel('E-W speed [m/s]'); datetick('x'); axis tight;

subplot(2,1,2); hold on

for i_b=1:numel(data.block.date)-1

fill([data.block.date(i_b) data.block.date(i_b) data.block.date(i_b+1) data.block.date(i_b+1)],[25 -25 -25 25],data.block.col(i_b,:))

end

fill([min(data.day) min(data.day) data.block.date(1) data.block.date(1)],[25 -25 -25 25],data.block.col(4,:))

fill([data.block.date(end) data.block.date(end) max(data.day) max(data.day) ],[25 -25 -25 25],data.block.col(4,:))

plot(fmulti.T(:),fmulti.Mvs(:),'.k');

for i_b=1:numel(data.block.date)

id = find(~fmulti.isnan & fmulti.B==i_b);

t_tmp = sort(unique(fmulti.T(id)));

plot(t_tmp, fmulti.f_trend(t_tmp,mean(data.lat(fmulti.R(id)))*ones(numel(t_tmp),1),mean(data.lon(fmulti.R(id)))*ones(numel(t_tmp),1), fmulti.beta_v(:,i_b)),'-r','Linewidth',2)

end

ylabel('N-S speed [m/s]'); datetick('x'); axis tight;

figure('position',[0 0 1000 400]); hold on;

[LAT, LON] = meshgrid(floor(min(data.lat)):.1:ceil(max(data.lat)),floor(min(data.lon)):.1:ceil(max(data.lon)));

for i_b=1:numel(data.block.date)

id = find(~fmulti.isnan & fmulti.B==i_b);

subplot(2,4,i_b,'Color',data.block.col(i_b,:)); hold on; h=worldmap([floor(min(data.lat)) ceil(max(data.lat))], [floor(min(data.lon)) ceil(max(data.lon))]);

setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')

surfm(LAT,LON,fmulti.f_trend(mean(fmulti.T(id))*ones(size(LAT)),LAT,LON,fmulti.beta_u(:,i_b)))

plotm(coastlat, coastlon,'k'); geoshow('worldrivers.shp','Color', 'blue')

[G,tmp]=findgroups(fmulti.R(id));

tmp2 = splitapply(@mean,fmulti.Mus(id),G);

scatterm(data.lat(tmp),data.lon(tmp),80,tmp2,'filled','MarkerEdgeColor','k');

title(data.block.date_str{i_b},'Color',data.block.col(i_b,:) ); caxis([min(tmp2) max(tmp2)]); colorbar

subplot(2,4,4+i_b,'Color',data.block.col(i_b,:)); hold on; h=worldmap([floor(min(data.lat)) ceil(max(data.lat))], [floor(min(data.lon)) ceil(max(data.lon))]);

setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')

surfm(LAT,LON,fmulti.f_trend(mean(fmulti.T(id))*ones(size(LAT)),LAT,LON,fmulti.beta_v(:,i_b)))

plotm(coastlat, coastlon,'k'); geoshow('worldrivers.shp','Color', 'blue')

tmp2 = splitapply(@mean,fmulti.Mvs(id),G);

scatterm(data.lat(tmp),data.lon(tmp),80,tmp2,'filled','MarkerEdgeColor','k');

title(data.block.date_str{i_b},'Color',data.block.col(i_b,:) ); caxis([min(tmp2) max(tmp2)]); colorbar

end

figure('position',[0 0 1000 400]);hold on;

[LAT, LON] = meshgrid(floor(min(data.lat)):1:ceil(max(data.lat)),floor(min(data.lon)):1:ceil(max(data.lon)));

for i_b=1:numel(data.block.date)

id = find(~fmulti.isnan & fmulti.B==i_b);

subplot(1,4,i_b,'Color',data.block.col(i_b,:)); hold on; h=worldmap([floor(min(data.lat)) ceil(max(data.lat))], [floor(min(data.lon)) ceil(max(data.lon))]);

setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')

plotm(coastlat, coastlon,'k'); geoshow('worldrivers.shp','Color', 'blue')

[G,tmp]=findgroups(fmulti.R(id));

tmpu = splitapply(@mean,fmulti.Mus(id),G);

tmpv = splitapply(@mean,fmulti.Mvs(id),G);

quiverm(LAT,LON,fmulti.f_trend(mean(fmulti.T(id))*ones(size(LAT)),LAT,LON,fmulti.beta_u(:,i_b)),fmulti.f_trend(mean(fmulti.T(id))*ones(size(LAT)),LAT,LON,fmulti.beta_v(:,i_b)),'k');

quiverm(data.lat(tmp),data.lon(tmp),tmpu,tmpv,'r');

title(data.block.date_str{i_b},'Color',data.block.col(i_b,:) ); caxis([min(tmp2) max(tmp2)]); colorbar

end

clear LAT LON T t_tmp tmp* G X id i_b h

Covariance Multi

cov.d=[0 1 150 200 300 400 500 750 1000 2000];
cov.t=[0 1 2 3 4 5 10 20];
[cov.D,cov.T] = meshgrid(cov.d(1:end-1)+diff(cov.d)/2,cov.t(1:end-1));
cov.D_emp=nan(size(cov.D,1),size(cov.D,2),numel(data.block.date)-1);
cov.T_emp = cov.D_emp;
cov.emp_grid_u = nan(size(cov.D));
cov.emp_grid_v = cov.emp_grid_u;
cov.nb=nan(numel(data.block.date)-1,1);
for i_b=1:numel(data.block.date)
    id = find(~fmulti.isnan & fmulti.B==i_b);
    cov.nb(i_b)=sum(id(:));
    Dtime=squareform(pdist(fmulti.T(id)));
    Ddist=squareform(pdist([data.lat(fmulti.R(id)) data.lon(fmulti.R(id))],@lldistkm));
    
    cov_emp_u = bsxfun(@times, fmulti.Mus_n(id)', fmulti.Mus_n(id));
    cov_emp_v = bsxfun(@times, fmulti.Mvs_n(id)', fmulti.Mvs_n(id));
    
    for i_t=1:numel(cov.t)-1
        for i_d=1:numel(cov.d)-1
            id = Dtime>=cov.t(i_t) & Dtime<cov.t(i_t+1) & Ddist>=cov.d(i_d) & Ddist<cov.d(i_d+1);
            cov.D_emp(i_t,i_d,i_b) = nanmean(Ddist(id));
            cov.T_emp(i_t,i_d,i_b) = nanmean(Dtime(id));
            cov.emp_grid_u(i_t,i_d,i_b) = nanmean( cov_emp_u(id) );
            cov.emp_grid_v(i_t,i_d,i_b) = nanmean( cov_emp_v(id) );
        end
    end
end

Figure

figure('position',[0 0 1000 400]);

for i_b = 1:numel(data.block.date)

subplot(2,4,i_b); surf(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),cov.emp_grid_u(:,:,i_b));

xlabel('Distance [km]'); ylabel('Time [days]'); title(data.block.date_str{i_b},'Color',data.block.col(i_b,:));

subplot(2,4,i_b,'color',data.block.col(i_b,:));

subplot(2,4,4+i_b); surf(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),cov.emp_grid_v(:,:,i_b));

xlabel('Distance [km]'); ylabel('Time [days]'); title(data.block.date_str{i_b},'Color',data.block.col(i_b,:));

subplot(2,4,4+i_b,'color',data.block.col(i_b,:));

end

Covariance function fit

Gneiting_fit = @(parm,i_b) parm(1)*(cov.D_emp(:,:,i_b)==0 & cov.T_emp(:,:,i_b)==0) + parm(2).*Gneiting(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),parm(3),parm(4),parm(5),parm(6),parm(7));
parm0 =[8 10 500 5 .5 .5 .5 ]';
for i_b=1:numel(data.block.date)
    id = find(~fmulti.isnan & fmulti.B==i_b);
    
    nugget_max = abs(diff(maxk(reshape(cov.emp_grid_u(:,:,i_b),1,[]),2)));
    parm_min=[0.8*nugget_max 1 100 1 .8 .8 0]';
    parm_max=[1.2*nugget_max nancov(fmulti.Mus_n(id)) 3000 10 1 1 1 ]';
    rmse = @(parm) nansum(nansum( ( 1000*max(cov.emp_grid_u(:,:,i_b),0.1).*(Gneiting_fit(parm,i_b) - cov.emp_grid_u(:,:,i_b)) ).^2 )); % 1000*max(cov.emp_grid(:,:,i_b),0).*
    fmulti.cov.parm_u(:,i_b) = fmincon(rmse,parm0,[],[],[1 1 0 0 0 0 0],nancov(fmulti.Mus_n(id)),parm_min,parm_max,[],optimoptions('fmincon','Display','off'));
    
    nugget_max = abs(diff(maxk(reshape(cov.emp_grid_v(:,:,i_b),1,[]),2)));
    parm_min=[0.9*nugget_max 1 100 1 .8 .8 0 ]';
    parm_max=[1.2*nugget_max nancov(fmulti.Mvs_n(id)) 3000 10 1 1 1 ]';
    rmse = @(parm) nansum(nansum( ( 1000*max(cov.emp_grid_v(:,:,i_b),0.1).*(Gneiting_fit(parm,i_b) - cov.emp_grid_v(:,:,i_b)) ).^2 )); % 1000*max(cov.emp_grid(:,:,i_b),0.1).*
    fmulti.cov.parm_v(:,i_b) = fmincon(rmse,parm0,[],[],[1 1 0 0 0 0 0],nancov(fmulti.Mvs_n(id)),parm_min,parm_max,[],optimoptions('fmincon','Display','off'));
end

fmulti.cov.C_un = @(D,T,parm) parm(1)*(D==0 & T==0) + parm(2).*Gneiting(D,T,parm(3),parm(4),parm(5),parm(6),parm(7));
fmulti.cov.C_vn = @(D,T,parm) parm(1)*(D==0 & T==0) + parm(2).*Gneiting(D,T,parm(3),parm(4),parm(5),parm(6),parm(7));
fmulti.cov.C_u  = @(D,T,parm) parm(2).*Gneiting(D,T,parm(3),parm(4),parm(5),parm(6),parm(7));
fmulti.cov.C_v  = @(D,T,parm) parm(2).*Gneiting(D,T,parm(3),parm(4),parm(5),parm(6),parm(7));

Ilustration of the fitted covariance function

[tmpD,tmpT] = meshgrid(cov.d(1):1:cov.d(end),cov.t(1):0.1:cov.t(end));

figure('position',[0 0 1000 400]);

for i_b = 1:numel(data.block.date)

subplot(2,4,i_b); hold on

s=surf(tmpD,tmpT,fmulti.cov.C_un(tmpD, tmpT,fmulti.cov.parm_u(:,i_b))); s.EdgeColor='none';

surf(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),cov.emp_grid_u(:,:,i_b));

xlabel('Distance [km]'); ylabel('Time [days]'); title(data.block.date_str{i_b},'Color',data.block.col(i_b,:));

plot3([0 0],[0 0],[fmulti.cov.parm_u(2,i_b) sum(fmulti.cov.parm_u(1:2,i_b))],'r', 'linewidth',2)

subplot(2,4,i_b,'color',data.block.col(i_b,:)); view(3)

subplot(2,4,4+i_b); hold on

s=surf(tmpD,tmpT,fmulti.cov.C_vn(tmpD, tmpT,fmulti.cov.parm_v(:,i_b))); s.EdgeColor='none';

surf(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),cov.emp_grid_v(:,:,i_b));

xlabel('Distance [km]'); ylabel('Time [days]'); title(data.block.date_str{i_b},'Color',data.block.col(i_b,:));

plot3([0 0],[0 0],[fmulti.cov.parm_v(2,i_b) sum(fmulti.cov.parm_v(1:2,i_b))],'r', 'linewidth',2)

subplot(2,4,4+i_b,'color',data.block.col(i_b,:)); view(3)

end

Clear all un-necessary variable

%clear C cov* Ddist Dtime Gneiting_fit h i* L* parm* rmse s tmp*

Intra-night scale

At the intra-night scale, we model the the fluctuation around the multi-night scale

fintra.Ius = data.us-fmulti.Mus(data.day_id,:);
fintra.Ivs = data.vs-fmulti.Mvs(data.day_id,:);
fintra_isnan = isnan(fintra.Ius);

Block decomposition

fintra_b = nan(size(data.time));
for i_b=1:numel(data.block.date)
    if i_b==numel(data.block.date)
        fintra_b(datenum(data.time)>=data.block.date(4) | datenum(data.time)<data.block.date(1)) = i_b;
    else
        fintra_b(datenum(data.time)>=data.block.date(i_b) & datenum(data.time)<data.block.date(i_b+1)) = i_b;
    end
end
fintra_B = repmat(fintra_b,1,data.nrad);

In addition to present a trend of the mean, the intra-night components also contains a non-stationarity of variance along the night (see figure below). We fit another polynomial function to normalize the variance

fintra_poly_degree=4;
X=bsxfun(@power,data.NNT(:),fintra_poly_degree:-1:0);
fintra.f_trend = @(t,beta_v) reshape(bsxfun(@power,t(:),fintra_poly_degree:-1:0) * beta_v(:),size(t));

We will then compute the detrended value of the multi-night scale componenent

fintra_Ius_var = nan(size(fintra.Ius));
fintra_Ivs_var = fintra_Ius_var;
fintra_beta_u = nan(fintra_poly_degree+1,numel(data.block.date));
fintra_beta_v = nan(fintra_poly_degree+1,numel(data.block.date));

We loop over all block...

for i_b=1:numel(data.block.date)
    id = fintra_B==i_b & ~fintra_isnan;
    fintra_beta_u(:,i_b) = (X(id,:)'*X(id,:))\(X(id,:)'*fintra.Ius(id).^2);
    fintra_beta_v(:,i_b) = (X(id,:)'*X(id,:))\(X(id,:)'*fintra.Ivs(id).^2);
    fintra_Ius_var(id) = fintra.f_trend(data.NNT(id),fintra_beta_u(:,i_b));
    fintra_Ivs_var(id) = fintra.f_trend(data.NNT(id),fintra_beta_v(:,i_b));
end

Illustration

NNT_axis=-1:.01:1;

figure('position',[0 0 1000 400]);

for i_b = 1:numel(data.block.date)

id = fintra_B==i_b & ~fintra_isnan;

subplot(2,4,i_b,'color',data.block.col(i_b,:)); hold on

plot(data.NNT(id),fintra.Ius(id).^2,'.k');

plot(NNT_axis, fintra.f_trend(NNT_axis,fintra_beta_u(:,i_b)),'-r','linewidth',2);

title(data.block.date_str{i_b},'color',data.block.col(i_b,:));

subplot(2,4,i_b,'color',data.block.col(i_b,:));

xlabel('Normalized Night Time (NNT)'); ylabel('var(B^p-M-\mu-\iota)');ylim([0 30]); xlim([-1 1])

subplot(2,4,4+i_b,'color',data.block.col(i_b,:)); hold on

plot(data.NNT(id),fintra.Ivs(id).^2,'.k');

plot(NNT_axis, fintra.f_trend(NNT_axis,fintra_beta_v(:,i_b)),'-r','linewidth',2);

title(data.block.date_str{i_b},'color',data.block.col(i_b,:));

subplot(2,4,4+i_b,'color',data.block.col(i_b,:));

xlabel('Normalized Night Time (NNT)'); ylabel('var(B^p-M-\mu-\iota)');ylim([0 30]); xlim([-1 1])

end

Finally, combining the polynomial of the curve and the correction for variance, we find the de-trended intra-night scale I

fintra.Ius_n = fintra.Ius(~fintra_isnan) ./ sqrt(fintra_Ius_var(~fintra_isnan));
fintra.Ivs_n = fintra.Ivs(~fintra_isnan) ./ sqrt(fintra_Ivs_var(~fintra_isnan));
clear X fintra_* NNT_axis i*

4.2 Covariance function of

Built the matrix of distances of the data. Distence (or difference) in lattitude, longitude, time and value.

Ddist_sm = squareform(pdist([data.lat data.lon], @lldistkm));
tmp=repmat(1:data.nrad,data.ntime,1);
fintra_radar = uint8(tmp(~fintra_isnan));
tmp=repmat(data.day(data.day_id),1,data.nrad);
fintra_day = single(tmp(~fintra_isnan));
tmp=repmat(datenum(data.time),1,data.nrad);
fintra_time = tmp(~fintra_isnan);
cov.d=[0 1 100 150 200 250 300 600 1000];
cov.t=[0 .01 .02 .05 .07 .1 .2];
[cov.D,cov.T]=meshgrid(cov.d(1:end-1)+diff(cov.d)/2,cov.t(1:end-1)+diff(cov.t)/2);
cov.D_emp=nan(size(cov.D,1),size(cov.D,2),numel(data.block.date)); 
cov.T_emp=nan(size(cov.D_emp));
cov.emp_grid_u=nan(size(cov.D_emp));
cov.emp_grid_v=nan(size(cov.D_emp));
for i_b=1:numel(data.block.date)
    D = nan(1000000000,2);
    Did = zeros(1000000000,2,'uint32');
    i=0;
    for i_t=data.day(data.day_b==i_b)'
        id = find(fintra_day==i_t);
        tmp1 = [reshape(Ddist_sm(fintra_radar(id),fintra_radar(id)),[],1) reshape(squareform(pdist(fintra_time(id))),[],1)];
        tmp2 = [uint32(repmat(id,numel(id),1)) uint32(repelem(id,numel(id)))];
        
        B = tril(ones(numel(id)));
        id2=B(:)==1 & tmp1(:,1)<cov.d(end) & tmp1(:,2)<cov.t(end);
        D(i+(1:sum(id2)),:)=tmp1(id2,:);
        Did(i+(1:sum(id2)),:)=tmp2(id2,:);
        i=i+sum(id2);
    end
    D=D(1:i,:);
    Did=Did(1:i,:);
    
    for i_d=1:numel(cov.d)-1
        id1 = find(D(:,1)>=cov.d(i_d) & D(:,1)<cov.d(i_d+1));
        for i_t=1:numel(cov.t)-1
            id2 = id1(D(id1,2)>=cov.t(i_t) & D(id1,2)<cov.t(i_t+1));
            cov.emp_grid_u(i_t,i_d,i_b)= mean(fintra.Ius_n(Did(id2,1)) .* fintra.Ius_n(Did(id2,2)));
            cov.emp_grid_v(i_t,i_d,i_b)= mean(fintra.Ivs_n(Did(id2,1)) .* fintra.Ivs_n(Did(id2,2)));
            cov.D_emp(i_t,i_d,i_b) = nanmean(D(id2,1));
            cov.T_emp(i_t,i_d,i_b) = nanmean(D(id2,2));
        end
    end
end

Figure

figure('position',[0 0 1000 400]);

for i_b = 1:numel(data.block.date)

subplot(2,4,i_b); surf(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),cov.emp_grid_u(:,:,i_b));

xlabel('Distance [km]'); ylabel('Time [days]'); title(data.block.date_str{i_b},'Color',data.block.col(i_b,:));

subplot(2,4,i_b,'color',data.block.col(i_b,:));

subplot(2,4,4+i_b); surf(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),cov.emp_grid_v(:,:,i_b));

xlabel('Distance [km]'); ylabel('Time [days]'); title(data.block.date_str{i_b},'Color',data.block.col(i_b,:));

subplot(2,4,4+i_b,'color',data.block.col(i_b,:));

end

Gneiting_fit = @(parm,i_b) parm(1)*(cov.D_emp(:,:,i_b)==0 & cov.T_emp(:,:,i_b)==0) + parm(2)*(cov.D_emp(:,:,i_b)==0) + parm(3).*Gneiting(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),parm(4),parm(5),parm(6),parm(7), parm(8));
parm0 =[0 0.2 1 100 0.1 .5 0 .5]';
parm_min=[0 0.3 0 0.0001 0.1 0 0.8 0]';
for i_b=1:numel(data.block.date)
    parm_max = [cov.emp_grid_u(1,1,i_b) cov.emp_grid_u(1,1,i_b) cov.emp_grid_u(1,1,i_b) 3000 10 1 1 1 ]';
    rmse = @(parm) nansum(nansum( double(cov.D>0.5) .* 1000.*max(cov.emp_grid_u(:,:,i_b),0).*( (Gneiting_fit(parm,i_b) - cov.emp_grid_u(:,:,i_b)).^2 ) )); % 1000*max(cov.emp_grid(:,:,i_b),0).*
    fintra.cov.parm_u(:,i_b) = fmincon(rmse,parm0,[],[],[1 1 1 0 0 0 0 0],cov.emp_grid_u(1,1,i_b),parm_min,parm_max,[],optimoptions('fmincon','Display','off'));
    
    parm_max = [cov.emp_grid_v(1,1,i_b) cov.emp_grid_v(1,1,i_b) cov.emp_grid_u(1,1,i_b) 3000 10 1 1 1 ]';
    rmse = @(parm) nansum(nansum( 1000*max(cov.emp_grid_v(:,:,i_b),0).*( (Gneiting_fit(parm,i_b) - cov.emp_grid_v(:,:,i_b)).^2 ) )); % 
    fintra.cov.parm_v(:,i_b) = fmincon(rmse,parm0,[],[],[1 1 1 0 0 0 0 0],cov.emp_grid_v(1,1,i_b),parm_min,parm_max,[],optimoptions('fmincon','Display','off'));
end

fintra.cov.C_un = @(D,T,parm) parm(1)*(D==0 & T==0) + parm(2)*(D==0) + parm(3).*Gneiting(D,T,parm(4),parm(5),parm(6),parm(7),parm(8));
fintra.cov.C_vn = @(D,T,parm) parm(1)*(D==0 & T==0) + parm(2)*(D==0) + parm(3).*Gneiting(D,T,parm(4),parm(5),parm(6),parm(7),parm(8));
fintra.cov.C_u  = @(D,T,parm) parm(3).*Gneiting(D,T,parm(4),parm(5),parm(6),parm(7),parm(8));
fintra.cov.C_v  = @(D,T,parm) parm(3).*Gneiting(D,T,parm(4),parm(5),parm(6),parm(7),parm(8));

Ilustration of the fitted covariance function

[tmpD,tmpT] = meshgrid(cov.d(1):5:cov.d(end),cov.t(1):0.001:cov.t(end));

figure('position',[0 0 1000 400]);

for i_b = 1:numel(data.block.date)

subplot(2,4,i_b); hold on

s=surf(tmpD,tmpT,fintra.cov.C_un(tmpD, tmpT, fintra.cov.parm_u(:,i_b))); s.EdgeColor='none';

surf(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),cov.emp_grid_u(:,:,i_b));

xlabel('Distance [km]'); ylabel('Time [days]'); title(data.block.date_str{i_b},'Color',data.block.col(i_b,:));

plot3([0 0],[0 0],[fintra.cov.parm_u(3,i_b) sum(fintra.cov.parm_u(2:3,i_b))],'r', 'linewidth',2)

plot3([0 0],[0 0],[sum(fintra.cov.parm_u(2:3,i_b)) sum(fintra.cov.parm_u(1:3,i_b))],'g', 'linewidth',2)

subplot(2,4,i_b,'color',data.block.col(i_b,:)); view(3)

subplot(2,4,4+i_b); hold on

s=surf(tmpD,tmpT,fintra.cov.C_vn(tmpD, tmpT, fintra.cov.parm_v(:,i_b))); s.EdgeColor='none';

surf(cov.D_emp(:,:,i_b),cov.T_emp(:,:,i_b),cov.emp_grid_v(:,:,i_b));

xlabel('Distance [km]'); ylabel('Time [days]'); title(data.block.date_str{i_b},'Color',data.block.col(i_b,:));

plot3([0 0],[0 0],[fintra.cov.parm_v(3,i_b) sum(fintra.cov.parm_v(2:3,i_b))],'r', 'linewidth',2)

plot3([0 0],[0 0],[sum(fintra.cov.parm_v(2:3,i_b)) sum(fintra.cov.parm_v(1:3,i_b))],'g', 'linewidth',2)

subplot(2,4,4+i_b,'color',data.block.col(i_b,:)); view(3)

end

5. Save

Save the parameters of the model.

save('data/Flight_inference.mat','fmulti','fintra','-v7.3')
% load('data/Flight_inference.mat')