Map Estimation

1. Load data

clear all; warning off
load('./data/Density_inference.mat','data');
load('./data/Density_estimationMap','g');
load('./data/Flight_inference.mat');
addpath('./functions/'); 
load coastlines.mat

4. Multi-night scale

Covariance matrix of weather radars

gfmulti.Mu_est = nan(g.nlm,g.nat);
gfmulti.Mu_sig = nan(g.nlm,g.nat);
gfmulti.Mv_est = nan(g.nlm,g.nat);
gfmulti.Mv_est = nan(g.nlm,g.nat);
for i_b=1:numel(data.block.date)
    id = ~fmulti.isnan & fmulti.B==i_b;
    Dtime=squareform(pdist(fmulti.T(id)));
    Ddist_aa=squareform(pdist([data.lat(fmulti.R(id)) data.lon(fmulti.R(id))],@lldistkm));
    
    Ddist_ab = pdist2([data.lat data.lon], [g.lat2D(g.latlonmask) g.lon2D(g.latlonmask)], @lldistkm);
    Ddist_sf = Ddist_ab(fmulti.R(id), repmat((1:g.nlm)',sum(g.day_b==i_b),1));
    Dtime_sf = pdist2(fmulti.T(id), repelem(datenum(g.day(g.day_b==i_b)),g.nlm,1));
    
    % u
    Caa = fmulti.cov.C_u(Ddist_aa,Dtime,i_b);
    Cab = fmulti.cov.C_u(Ddist_sf,Dtime_sf,i_b);
    Lambda = inv(Caa)*Cab;
    gfmulti.Mu_est(:,g.day_b==i_b) = reshape((Lambda' * fmulti.Mus(id)) + fmulti.f_trend(repelem(datenum(g.day(g.day_b==i_b)),g.nlm,1), repmat(g.lat2D(g.latlonmask),sum(g.day_b==i_b),1), repmat(g.lon2D(g.latlonmask),sum(g.day_b==i_b),1), fmulti.beta_u(:,i_b)), g.nlm, sum(g.day_b==i_b));
    gfmulti.Mu_sig(:,g.day_b==i_b) = reshape( sqrt(sum(fmulti.cov.parm_u(1:2,i_b)) - sum(Lambda.*Cab)), g.nlm, sum(g.day_b==i_b));
    
    % v
    Caa = fmulti.cov.C_v(Ddist_aa,Dtime,i_b);
    Cab = fmulti.cov.C_v(Ddist_sf,Dtime_sf,i_b);
    Lambda = inv(Caa)*Cab;
    gfmulti.Mv_est(:,g.day_b==i_b) = reshape((Lambda' * fmulti.Mvs(id)) + fmulti.f_trend(repelem(datenum(g.day(g.day_b==i_b)),g.nlm,1), repmat(g.lat2D(g.latlonmask),sum(g.day_b==i_b),1), repmat(g.lon2D(g.latlonmask),sum(g.day_b==i_b),1), fmulti.beta_v(:,i_b)), g.nlm, sum(g.day_b==i_b));
    gfmulti.Mv_sig(:,g.day_b==i_b) = reshape( sqrt(sum(fmulti.cov.parm_v(1:2,i_b)) - sum(Lambda.*Cab)), g.nlm, sum(g.day_b==i_b));
        
end

Figure

% tmp=nan(g.nlat, g.nlon);
% figure('position',[0 0 1000 600]); 
% i=1;
% for i_t = 1:12:4*6
%     subplot(4,6,i); hold on; i=i+1;
%     h=worldmap([min([g.lat]) max([g.lat])], [min([g.lon]) max([g.lon])]); 
%     setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
%     geoshow('landareas.shp', 'FaceColor', [215 215 215]./255); geoshow('worldrivers.shp','Color', 'blue')
%     tmp(g.latlonmask) = gfmulti.Mu_est(:,i_t);
%     surfm(g.lat2D,g.lon2D,tmp); caxis([-15 15])
%     i_tt=find(datenum(g.day(i_t))==data.day);
%     scatterm(data.lat(~isnan(fmulti.Mus(i_tt,:))),data.lon(~isnan(fmulti.Mus(i_tt,:))),100,fmulti.Mus(i_tt,~isnan(fmulti.Mus(i_tt,:))),'filled','MarkerEdgeColor','k')
%     title(datestr(g.day(i_t),'dd-mmm'))
% end

Export as gif

% % filename='data/Density_estimationMap_M';
% filename='data/Density_estimationMap_M';
% Frame(g.nat) = struct('cdata',[],'colormap',[]); 
% fig=figure(2); 
% h=worldmap([min([g.lat]) max([g.lat])], [min([g.lon]) max([g.lon])]); 
% setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
% geoshow('landareas.shp', 'FaceColor', [215 215 215]./255); geoshow('worldrivers.shp','Color', 'blue')
% set(gcf,'color','w'); hsurf=[];hscat=[];
% caxis([-3 3]); 
% tmp=nan(g.nlat, g.nlon);
% %caxis([0 80]);  c=colorbar;c.Label.String = 'Bird density [bird/km^2]';c.Label.FontSize=14;
% for i_t = 1:g.nat
%     delete(hsurf); delete(hscat)
%     tmp(g.latlonmask) = gmulti.M_est(:,i_t);
%     hsurf=surfm(g.lat2D,g.lon2D,tmp);
%     i_tt=find(datenum(g.day(i_t))==data.day);
%     hscat=scatterm(data.lat(~isnan(multi.M_m(i_tt,:))),data.lon(~isnan(multi.M_m(i_tt,:))),100,multi.M_m(i_tt,~isnan(multi.M_m(i_tt,:))),'filled','MarkerEdgeColor','k');
%     
%     title(datestr(g.day(i_t),'dd-mmm'),'FontSize',14); drawnow;
%     Frame(i_t) = getframe(fig);
%     [imind,cm] = rgb2ind(frame2im(Frame(i_t)),256); 
%     if i_t == 1
%         imwrite(imind,cm,[filename '.gif'],'gif','Loopcount',inf);
%     else
%         imwrite(imind,cm,[filename '.gif'],'gif','WriteMode','append');
%     end
% end

5. Intra-night scale

Covariance matrix of weather radar data

gfintra.Ius_est = nan(g.nlm,g.nt);
gfintra.Ius_sig = nan(g.nlm,g.nt);
gfintra.Ivs_est = nan(g.nlm,g.nt);
gfintra.Ivs_sig = nan(g.nlm,g.nt);
Ddist_rr_sm = squareform(pdist([data.lat data.lon], @lldistkm));
Ddist_rg_sm = pdist2([data.lat data.lon], [g.lat2D(g.latlonmask) g.lon2D(g.latlonmask)], @lldistkm);
tic
for i_t=1:g.nt
    disp([num2str(round(i_t/g.nat*100)) '%  ' num2str(toc/60) 'min'])
    i_b = g.day_b(i_t);
    id_r = fintra.day==datenum(g.day(i_t));
    id_g = g.day_id == i_t;
    
    Dtime_rr = squareform(pdist(fintra.time(id_r)));
    Ddist_rr = Ddist_rr_sm(fintra.radar(id_r),fintra.radar(id_r));
    
    Dtime_rg = pdist2(fintra.time(id_r),datenum(g.time(id_g)));
    Dtime_rg = Dtime_rg(:,repelem((1:sum(id_g))',g.nlm,1));
    Ddist_rg = Ddist_rg_sm(fintra.radar(id_r),:);
    Ddist_rg = repmat(Ddist_rg,1,sum(id_g));
    
    id = g.mask_day(:,id_g);
    id2 = false(size(gfintra.Ius_est));
    id2(:,id_g)=id;
    
    % u
    Crr = fintra.cov.C_u(Ddist_rr,Dtime_rr,i_b);
    Crg = fintra.cov.C_u(Ddist_rg(:,id(:)),Dtime_rg(:,id(:)),i_b);
    Lambda = inv(Crr)*Crg;
    gfintra.Ius_est(id2) = (Lambda' * fintra.Ius_n(id_r)) .* sqrt(fintra.f_trend(g.NNT(id2),fintra.beta_u(:,i_b)));
    gfintra.Ius_sig(id2) = sqrt( (sum(fintra.cov.parm_u(1:2,i_b)) - sum(Lambda.*Crg))  .* fintra.f_trend(g.NNT(id2)',fintra.beta_u(:,i_b)) ); 
    
    % v
    Crr = fintra.cov.C_v(Ddist_rr,Dtime_rr,i_b);
    Crg = fintra.cov.C_v(Ddist_rg(:,id(:)),Dtime_rg(:,id(:)),i_b);
    Lambda = inv(Crr)*Crg;
    gfintra.Ivs_est(id2) = (Lambda' * fintra.Ivs_n(id_r)) .* sqrt(fintra.f_trend(g.NNT(id2),fintra.beta_v(:,i_b)));
    gfintra.Ivs_sig(id2) = sqrt( (sum(fintra.cov.parm_v(1:2,i_b)) - sum(Lambda.*Crg))  .* fintra.f_trend(g.NNT(id2)',fintra.beta_v(:,i_b)) ); 
        
end
63%  0.00020022min
63%  1.3447min
63%  2.6335min
64%  4.1346min
64%  5.3797min
64%  6.4421min
65%  7.9752min
65%  9.595min
65%  11.3633min
66%  12.9379min
66%  14.4416min
66%  15.9009min
66%  17.346min
67%  18.2293min
67%  19.1239min
67%  20.0114min
68%  20.8113min
68%  21.6486min
68%  22.5732min
69%  24.5013min
69%  25.6634min
69%  26.2404min
70%  27.3124min
70%  29.3997min
70%  31.2628min
70%  33.1866min
71%  35.2206min
71%  37.37min
71%  39.4581min
72%  41.2183min
72%  42.4299min
72%  43.8838min
73%  46.1516min
73%  48.2846min
73%  50.9065min
74%  52.0577min
74%  54.4468min
74%  57.0572min
75%  59.7477min
75%  61.8377min
75%  63.4283min
75%  65.3883min
76%  67.0533min
76%  68.431min
76%  70.5588min
77%  73.3181min
77%  75.873min
77%  78.5863min
78%  81.7525min
78%  84.5773min
78%  87.3403min
79%  89.9222min
79%  91.4315min
79%  93.1336min
80%  95.9779min
80%  96.9507min
80%  99.0176min
80%  100.4399min
81%  101.5479min
81%  103.3783min
81%  105.2804min
82%  107.1296min
82%  110.2455min
82%  113.6453min
83%  117.1945min
83%  119.5047min
83%  121.8069min
84%  123.8302min
84%  125.5872min
84%  127.1989min
84%  127.8508min
85%  129.3748min
85%  130.9591min
85%  134.5831min
86%  138.1884min
86%  142.579min
86%  145.9066min
87%  148.4292min
87%  150.4882min
87%  152.0014min
88%  153.4158min
88%  156.0358min
88%  158.3165min
89%  159.5576min
89%  161.0168min
89%  162.4561min
89%  164.0192min
90%  165.0165min
90%  166.0615min
90%  167.2492min
91%  169.1244min
91%  169.4264min
91%  169.6019min
92%  170.4447min
92%  172.843min
92%  173.405min
93%  174.3812min
93%  174.9055min
93%  175.0369min
93%  176.0008min
94%  177.9381min
94%  180.224min
94%  182.7225min
95%  184.9022min
95%  186.9295min
95%  187.7318min
96%  188.7592min
96%  190.9804min
96%  191.8958min
97%  192.5919min
97%  193.0609min
97%  193.7238min
98%  194.887min
98%  195.9352min
98%  197.4169min
98%  199.8948min
99%  202.1372min
99%  204.1253min
99%  205.9287min
100%  207.321min
100%  208.4891min

6 Reassemble

guv.u_est = gfmulti.Mu_est(:,g.day_id) + gfintra.Ius_est;
guv.v_est = gfmulti.Mv_est(:,g.day_id) + gfintra.Ivs_est;
guv.u_sig = sqrt(gfmulti.Mu_sig(:,g.day_id).^2 + gfintra.Ius_sig.^2);
guv.v_sig = sqrt(gfmulti.Mv_sig(:,g.day_id).^2 + gfintra.Ivs_sig.^2);

Figure

% fig=figure(2);  
% filename='data/Flight_estimationMap';
% Frame(g.nt) = struct('cdata',[],'colormap',[]);
% h=worldmap([min([g.lat]) max([g.lat])], [min([g.lon]) max([g.lon])]); 
% setm(h,'grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
% geoshow('landareas.shp', 'FaceColor', [215 215 215]./255); geoshow('worldrivers.shp','Color', 'blue')
% set(gcf,'color','w'); hsurf=[];hscat=[];
% caxis([-15 15]);
% tmp_u=nan(g.nlat, g.nlon); tmp_v=nan(g.nlat, g.nlon);
% for i_t = 1:g.nt
%     delete(hsurf); delete(hscat);
%     tmp_u(g.latlonmask)=guv.u_est(:,i_t);
%     tmp_v(g.latlonmask)=guv.v_est(:,i_t);
%     hsurf=quiverm(g.lat2D,g.lon2D,tmp_u,tmp_v);
% 
%     % id = abs(intra.time-datenum(g.time(i_t)))<1/24/4;
%     % hscat=scatterm(data.lat(intra.radar(id)),data.lon(intra.radar(id)),[],intra.I_m(id),'filled','MarkerEdgeColor','k');
% %     id = abs(data.time-g.time(i_t))<1/24/4;
% %     hscat=scatterm(data.lat,data.lon,[],nanmean(data.us(id,:),1),'filled','MarkerEdgeColor','k');
%     
%     title(datestr(g.time(i_t),'dd-mmm HH:MM'),'FontSize',14); drawnow
% %     Frame(i_t) = getframe(fig);
% %     [imind,cm] = rgb2ind(frame2im(Frame(i_t)),256); 
% %     if i_t == 2
% %         imwrite(imind,cm,[filename '.gif'],'gif', 'DelayTime',0.1,'Loopcount',inf);
% %     else
% %         imwrite(imind,cm,[filename '.gif'],'gif', 'DelayTime',0.1,'WriteMode','append');
% %     end
% keyboard
% end
% v=VideoWriter([filename '.avi']);
% v.FrameRate = 4;
% v.Quality = 75;
% open(v);
% writeVideo(v, Frame);
% close(v);

7. Save

save('data/Flight_estimationMap','guv','-v7.3');
% load('data/Flight_estimationMap')