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