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')