Modeling Mird Migration as a Flow
Table of Contents
Modeling Mird Migration as a Flow
Methodology
Flow model
Smoothing
Migratory Processes
Results
Nightly migratory pulses and stopovers
Year-round accumulation of migratory birds on the ground
Seasonal Flow: In and Out
Seasonal Flow: Within the study area
Appendix
Illutration main text
Supplementary Material not published
Methodology
Flow model
Smoothing
Migratory Processes
Results
Nightly migratory pulses and stopovers
Year-round accumulation of migratory birds on the ground
Seasonal Flow: In and Out
Seasonal Flow: Within the study area
Appendix
Illutration main text
Supplementary Material not published
Methodology
Flow model
The equation of continuity is
Load data
clear all;
load('data/Density_inference.mat','data');
load('data/Density_estimationMap','g','gd');
load('data/Flight_estimationMap','guv');
addpath('./functions/')
load coastlines;
% col2=colormap(brewermap([],'Spectral'));
col2=colormap(brewermap([],'Paired'));
Build extended grid
Build a grid gext with a cell more on each side of the domain (ext=extanded) in order to compute the fluxes in/out of the existing grid
gext.lat = fillmissing([nan(1,1) ;g.lat ;nan(1,1)],'linear');
gext.lon = fillmissing([nan(1,1); g.lon; nan(1,1)],'linear');
gext.nlat = numel(gext.lat);
gext.nlon = numel(gext.lon);
[gext.lat2D, gext.lon2D] = ndgrid(gext.lat,gext.lon);
The extended grid also include the time between two time step of the previous grid. Yet, because there were no time step between the last time step of the end of night (sunrise) and the first of the next day (sunset), we add one at noon every day
gext.time = sort( [g.time; g.day-.5]);
Compute different index and other stuff usefull later for the extended grid
gext.nt = numel(gext.time);
gext.day_id=nan(gext.nt,1);
gext.day_id(ismember(gext.time,g.time)) = g.day_id;
gext.day_id(isnan(gext.day_id)) = gext.day_id(find(isnan(gext.day_id))+1);
tmp=nan(g.nlm,gext.nt);
tmp(:,ismember(gext.time,g.time)) = g.NNT;
gext.NNT=nan(gext.nlat,gext.nlon,gext.nt);
gext.NNT(padarray(repmat(g.latlonmask,1,1,gext.nt),[1 1 0], false))=tmp;
gext.mask_day = logical(gext.NNT>-1 & gext.NNT<1);
gext.time_b=nan(gext.nt,1);
gext.time_b(ismember(gext.time,g.time)) = g.time_b;
gext.latlonmask = padarray(g.latlonmask,[1 1],false);
Compute grid resolution
dy = lldistkm([g.lat(1) g.lon(1)],[g.lat(2) g.lon(1)]);
dx = lldistkm([g.lat2D(:,1) g.lon2D(:,1)],[g.lat2D(:,1) g.lon2D(:,2)]);
area = repmat(dx*dy,1,g.nlon);
dxext = fillmissing([NaN;dx;NaN],'linear');
areaext = fillmissing(padarray(area,[1 1],nan),'linear',1);
areaext = fillmissing(areaext,'linear',2);
The time resolution is 15 min -> 15/60 hr
dt=15/60;
max(guv.u_est(:)/1000*60*60)*dt/min(dx) + max(guv.v_est(:)/1000*60*60)*dt/dy
Compute flight speed and density
Reshape the flight speed and density on a 3D grid (it is initially saved on a 2D grid, using only location with data)
vy = nan(g.nlat, g.nlon, gext.nt); vx=vy;
tmp = repmat(g.latlonmask,1,1,gext.nt); tmp(:,:,ismember(gext.time,g.day-.5))=false;
vy(tmp) = guv.v_est/1000*60*60; % m/s -> km/h (+) north, (-) south
vx(tmp) = guv.u_est/1000*60*60; % m/s -> km/h (+) east, (-) west
tmp2 = gd.dens_est; % bird/km^2 .* repmat(area(g.latlonmask),1,g.nt); % bird
Removing data when rain.
% tmp2(g.rain(repmat(g.latlonmask,1,1,g.nt))>data.mask_rain_thr) = 0;
Rho is not density be Number of birds here.
rho = nan(g.nlat, g.nlon, gext.nt);
rho(tmp) = tmp2; % bird/km^2
We also create the same tho wihtout nan (able to comupte difference)
rho0 = rho;
rho0(isnan(rho)) = 0; % bird/km^2
Compute the flux 
at +/- 1/2 grid cell.
at +/- 1/2 grid cell.First, add a nan layer in lat or lon direction
Phiy_pad = padarray(rho .* vy .* repmat(dx,1,g.nlon,gext.nt) ,[1 0 0],nan); % bird/km^2 * km/h *km -> bird/h
Phix_pad = padarray(rho .* vx .* dy,[0 1 0],nan);
and then, compute the flux at +/- 1/2 even if either previous or next cells is nan.
Phiy_h = movmean(Phiy_pad,[0 1],1,'omitnan','Endpoints','discard');
Phix_h = movmean(Phix_pad,[0 1],2,'omitnan','Endpoints','discard');
Compute the delta flux / delta distance 
First, add 0 padding for the outer zone (allows to compute boundary cell)
Phiy_h_0=padarray(Phiy_h,[1 1 0],0);
Phix_h_0=padarray(Phix_h,[1 1 0],0);
Then, replace the nan by zero to compute the change of fluxes even no fluxes on the other side (i.e. boundary cell)
Phiy_h_0(isnan(Phiy_h_0))=0;
Phix_h_0(isnan(Phix_h_0))=0;
Finally, compute the delta flux over delta distance.
dPhiy = diff(Phiy_h_0,1,1);
dPhix = diff(Phix_h_0,1,2);
Compute the variation of bird in/out each cell
F = (dPhiy + dPhix ).*dt; % bird/h * hr -> bird
Crank-nicolson Move F at t+1/2:
% F(F==0)=nan;
% F = movmean( F ,[0 1],3,'omitnan');
% F(isnan(F))=0;
F comprise both the intral change of bird and the bird going out of the area
Inner flux: remove the boundary cells
Fin = F;
%Fin(~gext.mask_day)=0;
Fin(repmat(~gext.latlonmask,1,1,gext.nt)) = 0;
Outer Flux. Note that Philat_h_0 and Philat_h_0 are 0 for the cell outside the domain, such at dPhilatdlat is equal to the flux at the outer cell (with the sign corresponding to the direction).
Fout.Fout = F;
%Fout.Fout(gext.mask_day)=0;
Fout.Fout(repmat(gext.latlonmask,1,1,gext.nt)) = 0;
Compute the landing/derpature W
rho0 is used (NaN replaced by 0) to account for landing as sunset and sunrise. Unit of W is [bird =bird (*15min/15min)], when (+) Birds entering the airspace (takingoff) and (-) Birds leaving the raispace (landing)
W.W = diff(rho0 .* repmat(area,1,1,gext.nt),1,3) + Fin(2:end-1,2:end-1,1:end-1);
Smoothing
W: standard deviation of 1% 3x0.25° = 0.75° -> ~50-80km
% std_deg = 1;
% ntab = 15;
% B = fspecial('gaussian', [ntab*2+1 ntab*2+1], std_deg);
% %B = fspecial('disk', 3);
% WW0 = padarray(g.latlonmask,[ntab ntab],'post');
%
% fftB=fft2(circshift(padarray(B,size(WW0)-size(B),'post'),-floor(size(B)/2)));
% N=ifft2(fft2(WW0).*fftB,'symmetric');
% N2=ifft2(fft2(~WW0).*fftB,'symmetric');
% N(N==0)=100; % whatever value
% W.Ws = ifft2(fft2(padarray(W.W,[ntab ntab 0],'post')).*fftB,'symmetric')+N2.*padarray(W.W,[ntab ntab 0],'post');
% W.Ws = W.Ws(1:(end-ntab),1:(end-ntab),:);
% W.Ws(repmat(~g.latlonmask,1,1,gext.nt-1))=0;
F: standard deviation of 1% 3x0.25° = 0.75° -> ~50-80km
% std_deg = 5;
% ntab = 15;
% B = fspecial('gaussian', [ntab*2+1 ntab*2+1], std_deg);
% %B = fspecial('disk', 3);
% FoutFout0 = padarray(~all(Fout.Fout==0,3),[ntab ntab],'post');
%
% fftB=fft2(circshift(padarray(B,size(FoutFout0)-size(B),'post'),-floor(size(B)/2)));
% N=ifft2(fft2(FoutFout0).*fftB,'symmetric');
% N2=ifft2(fft2(~FoutFout0).*fftB,'symmetric');
% N(N==0)=100; % whatever value
% Fout.Fouts = ifft2(fft2(padarray(Fout.Fout,[ntab ntab 0],'post')).*fftB,'symmetric')+N2.*padarray(Fout.Fout,[ntab ntab 0],'post');
% Fout.Fouts = Fout.Fouts(1:(end-ntab),1:(end-ntab),:);
% Fout.Fouts(repmat(all(Fout.Fout==0,3),1,1,gext.nt-1))=0;
%
% figure;
% subplot(1,2,1);
% imagesc(Fout.Fout(:,:,300)); caxis([-300 300])
% subplot(1,2,2);
% imagesc(Fout.Fouts(:,:,300)); caxis([-300 300])
Check methodology
Check for a step
i=20; i=70; i=56;% i=1095;
i=2605;i=2694;
gext.time(i);
figure('position',[0 0 1200 600]);
subplot(2,3,1); imagesc(Phix_pad(:,:,i),'AlphaData',~isnan(Phix_pad(:,:,i))); set(gca,'ydir','normal'); colorbar; title('\Phi_{x}')
subplot(2,3,4); imagesc(Phiy_pad(:,:,i),'AlphaData',~isnan(Phiy_pad(:,:,i))); set(gca,'ydir','normal'); colorbar; title('\Phi_{y}')
subplot(2,3,2); imagesc(dPhix(:,:,i),'AlphaData',~(0==dPhix(:,:,i))); set(gca,'ydir','normal'); colorbar; title('\Delta\Phi_{x}/\Delta x')
subplot(2,3,5); imagesc(dPhiy(:,:,i),'AlphaData',~(0==dPhiy(:,:,i))); set(gca,'ydir','normal'); colorbar; title('\Delta\Phi_{y}/\Delta y')
subplot(2,3,6); imagesc(F(:,:,i),'AlphaData',~(0==F(:,:,i))); set(gca,'ydir','normal'); colorbar; title('F=(\Delta\Phi_{y}+\Delta\Phi_{x})\Delta t')
Check the spatial aggregated mass balence.
a = reshape(nansum(nansum(Fout.Fout,1),2),1,[]);
b = reshape(nansum(nansum(rho0 .* repmat(area,1,1,gext.nt),1),2),1,[]);
c = reshape(nansum(nansum(W.W,1),2),1,[]);
figure('position',[0 0 1000 200]); hold on; xlabel('time');ylabel('Error')
plot(b(2:end) - b(1:end-1) - c - a(1:end-1))
See a single timestep.
figure('position',[0 0 1300 600]);
subplot(2,3,1); imagesc( rho0(:,:,i),'AlphaData',~(0==rho0(:,:,i))); title('\rho(t)');set(gca,'ydir','normal'); colorbar;
subplot(2,3,2); imagesc( rho0(:,:,i+1),'AlphaData',~(0==rho0(:,:,i+1))); title('\rho(t+1)'); set(gca,'ydir','normal'); colorbar;
subplot(2,3,3); imagesc( rho0(:,:,i+1)-rho0(:,:,i),'AlphaData',~(0==rho0(:,:,i))); title('\rho(t+1) - \rho(t)'); set(gca,'ydir','normal'); colorbar;
subplot(2,3,4); imagesc(Fin(:,:,i),'AlphaData',~(0==Fin(:,:,i))); title('F_{in}'); set(gca,'ydir','normal'); colorbar;
subplot(2,3,5); imagesc(W.W(:,:,i),'AlphaData',~(0==W.W(:,:,i))); title('W^{t\rightarrow t+1}=\rho(t+1) - \rho(t)+F_{in}'); set(gca,'ydir','normal'); colorbar;
subplot(2,3,6); imagesc(Fout.Fout(:,:,i),'AlphaData',~(0==Fout.Fout(:,:,i))); title('F_{out}'); set(gca,'ydir','normal'); colorbar;
Migratory Processes
Seperature the landing/landing (-) and takingoff/take-off (+)
W.landing=W.W; W.takingoff=W.W;
W.landing(W.landing>=0)=0;
W.landing2 = W.landing;
W.landing2(~gext.mask_day(2:end-1,2:end-1,1:end-1))=0;
W.takingoff(W.takingoff<=0)=0;
Fout.leaving=Fout.Fout; Fout.entering=Fout.Fout;
Fout.leaving(Fout.leaving>=0)=0;
Fout.entering(Fout.entering<=0)=0;
Compute daily sum
W.day.takingoff=nan(g.nlat,g.nlon,g.nat);
W.day.landing=nan(g.nlat,g.nlon,g.nat);
W.day.landing2=nan(g.nlat,g.nlon,g.nat);
Fout.day.entering=nan(gext.nlat,gext.nlon,g.nat);
Fout.day.leaving=nan(gext.nlat,gext.nlon,g.nat);
F_day=nan(gext.nlat,gext.nlon,g.nat);
for i_t=1:g.nat-1
idt=gext.day_id(1:end-1)==i_t;
W.day.takingoff(:,:,i_t)=nansum(W.takingoff(:,:,idt),3);
W.day.landing(:,:,i_t)=nansum(W.landing(:,:,idt),3);
W.day.landing2(:,:,i_t)=nansum(W.landing2(:,:,idt),3);
Fout.day.entering(:,:,i_t)=nansum(Fout.entering(:,:,idt),3);
Fout.day.leaving(:,:,i_t)=nansum(Fout.leaving(:,:,idt),3);
F_day(:,:,i_t)=nansum(F(:,:,idt),3);
end
W.day.takingoff(W.day.takingoff==0)=nan;
W.day.landing(W.day.landing==0)=nan;
W.day.landing2(W.day.landing2==0)=nan;
W.day.W=W.day.takingoff+W.day.landing;
Only keep with sufficient data.
knownday = ismember(datenum(g.day), unique(data.day(data.day_id)));
W.day.takingoff(:,:,~knownday)=nan;
W.day.landing(:,:,~knownday)=nan;
W.day.W(:,:,~knownday)=nan;
Compute timeseries
Ts.W = reshape(nansum(nansum(W.W,1),2),1,[]);
Ts.landing = reshape(nansum(nansum(W.landing,1),2),1,[]);
Ts.takingoff = reshape(nansum(nansum(W.takingoff,1),2),1,[]);
Ts.Fout.entering = reshape(nansum(nansum(Fout.entering,1),2),1,[]);
Ts.Fout.leaving = reshape(nansum(nansum(Fout.leaving,1),2),1,[]);
Ts.day.landing = reshape(nansum(nansum(W.day.landing,1),2),1,[]);
Ts.day.takingoff = reshape(nansum(nansum(W.day.takingoff,1),2),1,[]);
Ts.Fout.day.entering = reshape(nansum(nansum(Fout.day.entering,1),2),1,[]);
Ts.Fout.day.leaving = reshape(nansum(nansum(Fout.day.leaving,1),2),1,[]);
Ts.day.landing(~knownday)=nan;
Ts.day.takingoff(~knownday)=nan;
Ts.Fout.day.leaving(~knownday)=nan;
Ts.Fout.day.entering(~knownday)=nan;
Save
save('data/SinkSource.mat','W','Fout','Ts','gext','vy','vx','rho');
%load('data/SinkSource.mat')
Remove unused data
clear ans B dP* fftB N N2 ntab Phi* tmp* WW0
Results
Nightly migratory pulses and stopovers
figure('position',[0 0 1000 600]);
i_t=52:55;
for i_tt=1:numel(i_t)
subplot(3,numel(i_t),i_tt); hold on
h = worldmap([g.lat(1) g.lat(end)], [g.lon(1) g.lon(end)]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
geoshow('landareas.shp', 'FaceColor', [77 77 77]./255); geoshow('worldrivers.shp','Color', 'blue')
W_tmp=nansum(W.day.takingoff(:,:,i_t(i_tt)),3);
W_tmp(W_tmp==0)=NaN;
surfm(g.lat2D,g.lon2D,W_tmp./area );
% plotm(nansum(W_tmp(:) .* g.lat2D(:)) ./ nansum(W_tmp(:)), nansum(W_tmp(:) .* g.lon2D(:)) ./ nansum(W_tmp(:)),'.k')
title(datestr(g.day(i_t(i_tt))))
colorbar
caxis([-80 80])
colormap(gca,brewermap([],'Spectral'))
subplot(3,numel(i_t),numel(i_t)+i_tt); hold on
h = worldmap([g.lat(1) g.lat(end)], [g.lon(1) g.lon(end)]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
geoshow('landareas.shp', 'FaceColor', [77 77 77]./255); geoshow('worldrivers.shp','Color', 'blue')
gd_tmp = nan(g.nlat, g.nlon);
gd_tmp(g.latlonmask)=nanmean(gd.dens_est(:,i_t(i_tt)==g.day_id),2);
gd_tmp(gd_tmp==0)=NaN;
surfm(g.lat2D,g.lon2D,gd_tmp );
% plotm(data.lat,data.lon,'.k')
gu_tmp = nan(g.nlat, g.nlon); gv_tmp = nan(g.nlat, g.nlon);
gu_tmp(g.latlonmask)=nanmean(guv.u_est(:,i_t(i_tt)==g.day_id),2);
gu_tmp(gu_tmp==0)=NaN;
gv_tmp(g.latlonmask)=nanmean(guv.v_est(:,i_t(i_tt)==g.day_id),2);
gv_tmp(gv_tmp==0)=NaN;
quiverm(g.lat2D(1:5:end,1:5:end),g.lon2D(1:5:end,1:5:end),gu_tmp(1:5:end,1:5:end),gv_tmp(1:5:end,1:5:end),'k')
colorbar; %colormap(gca,brewermap([],'Spectral'))
caxis([0 50])
subplot(3,numel(i_t),2*numel(i_t)+i_tt); hold on
h = worldmap([g.lat(1) g.lat(end)], [g.lon(1) g.lon(end)]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
geoshow('landareas.shp', 'FaceColor', [77 77 77]./255); geoshow('worldrivers.shp','Color', 'blue')
W_tmp=nansum(W.day.landing2(:,:,i_t(i_tt)),3);
W_tmp(W_tmp==0)=NaN;
surfm(g.lat2D,g.lon2D,W_tmp./area ); %plotm(data.lat,data.lon,'.k')
% plotm(nansum(W_tmp(:) .* g.lat2D(:)) ./ nansum(W_tmp(:)), nansum(W_tmp(:) .* g.lon2D(:)) ./ nansum(W_tmp(:)),'.k')
colorbar
caxis([-80 80])
colormap(gca,brewermap([],'Spectral'))
end
% set(gcf,'color',[77 77 77]./255)
Map-daily
fig=figure;
ax = axes('Parent',fig,'units', 'normalized','position',[.05 0.15 .9 .8]);
ax1 = subplot(1,2,1);
h = worldmap([g.lat(1) g.lat(end)], [g.lon(1) g.lon(end)]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
geoshow(ax1, 'landareas.shp', 'FaceColor', [215 215 215]./255); geoshow('worldrivers.shp','Color', 'blue')
h1=surfm(g.lat2D,g.lon2D, W.day.takingoff(:,:,20)); colorbar;caxis([0 50000])
plotm(data.lat,data.lon,'.k')
ax2 = subplot(1,2,2);
h = worldmap([g.lat(1) g.lat(end)], [g.lon(1) g.lon(end)]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
geoshow(ax2, 'landareas.shp', 'FaceColor', [215 215 215]./255); geoshow('worldrivers.shp','Color', 'blue')
h2=surfm(g.lat2D,g.lon2D, -W.day.landing(:,:,20)); colorbar; caxis([0 50000])
plotm(data.lat,data.lon,'.k')
b = uicontrol('Parent',fig,'Style','slider','units', 'normalized','Position',[0.05 0.05 .9 0.05],...
'value',1, 'min',1, 'max',size(W.day.W,3),'SliderStep', [1/size(W.day.W,3), 0.1]);
b.Callback = @(es,ed) cellfun( @(x) feval(x,es,ed), {...
@(es,ed) set(h1,'CData', W.day.takingoff(:,:,round(es.Value+1))), ...
@(es,ed) set(h2,'CData', -W.day.landing(:,:,round(es.Value))), ...
@(es,ed) colorbar,...
@(es,ed) title(ax1,['takingoff :' datestr(g.day(round(es.Value+1)))]),...
@(es,ed) title(ax2,['landing :' datestr(g.day(round(es.Value)))]),...
});
Year-round accumulation of migratory birds on the ground
Times serie daily
figure('position',[0 0 1000 1000]); subplot(3,1,[1 2]); box on; grid on
hold on; ylabel('Number of birdss (millions)')
% min_f = 250; max_f= 200;
% 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)],[max_f -min_f -min_f max_f],data.block.col(i_b,:),'FaceAlpha',0.5,'EdgeColor','none')
% end
% fill([min(data.day) min(data.day) data.block.date(1) data.block.date(1)],[max_f -min_f -min_f max_f],data.block.col(4,:),'FaceAlpha',0.5,'EdgeColor','none')
% fill([data.block.date(end) data.block.date(end) max(data.day) max(data.day) ],[max_f -min_f -min_f max_f],data.block.col(4,:),'FaceAlpha',0.5,'EdgeColor','none')
fill(datenum([g.day(1); g.day; g.day(end)]),[0 cumsum((Ts.Fout.day.leaving+Ts.Fout.day.entering)/1000000,'omitnan') 0],'k','FaceAlpha',0.2);
h1=plot(datenum(g.day),cumsum((Ts.Fout.day.leaving+Ts.Fout.day.entering)/1000000,'omitnan'),'-k','linewidth',2 ,'DisplayName','Cumulative \Delta takingoff/landing = # on the ground');
h2=bar(datenum(g.day),Ts.day.landing/1000000,'FaceColor',col2(6,:),'DisplayName','Landing');
h3=bar(datenum(g.day),Ts.day.takingoff/1000000,'FaceColor',col2(2,:),'DisplayName','Takingoff');
h4=bar(datenum(g.day),Ts.Fout.day.leaving/1000000,'FaceColor',col2(8,:),'DisplayName','Leaving');
h5=bar(datenum(g.day),Ts.Fout.day.entering/1000000,'FaceColor',col2(10,:),'DisplayName','Entering +');
h6=bar(datenum(g.day),(Ts.Fout.day.leaving+Ts.Fout.day.entering)/1000000,'k','DisplayName','\Delta takingoff-landing = \Delta Leaving- Entering');
datetick('x'); axis tight; legend([h2 h3 h4 h5 h6 h1]);
subplot(3,1,3); box on; grid on
hold on; ylabel('Number of birdss (millions)')
h1=fill(datenum([g.day ;flipud(g.day)]),-[cumsum((Ts.Fout.day.leaving)/1000000,'omitnan') fliplr(cumsum(-(Ts.Fout.day.entering)/1000000,'omitnan'))],'k','FaceAlpha',0.2,'DisplayName','Difference');
h2=plot(datenum(g.day),cumsum(-(Ts.Fout.day.leaving)/1000000,'omitnan'),'Color',col2(8,:),'linewidth',2 ,'DisplayName','Cumulative Leaving');
h3=plot(datenum(g.day),cumsum((Ts.Fout.day.entering)/1000000,'omitnan'),'Color',col2(10,:),'linewidth',2 ,'DisplayName','Cumulative Entering');
datetick('x'); xlabel('Date (2018)'); axis tight; legend([h2 h3 h1]);
Maximum Number of bird in a night
[max_takeoff,tmp] = max(Ts.day.takingoff(g.day<datetime('15-July-2018')));
disp(['Maximum number of bird taking-off is: ' num2str(max_takeoff/1000000) ' happening on the ' datestr(g.day(tmp)) '.'])
[max_takeoff,tmp] = max(Ts.day.takingoff);
disp(['Maximum number of bird taking-off is: ' num2str(max_takeoff/1000000) ' happening on the ' datestr(g.day(tmp)) '.'])
How many night to reach 50% of all derpature
tmp = Ts.day.takingoff(g.day<datetime('15-July-2018') & ~isnan(Ts.day.takingoff'));
disp(['50% of all bird taking-off in Spring happened in ' num2str(sum( (cumsum(sort(tmp,'descend'),'omitnan') ./ nansum(tmp)) < .5) )])
tmp = Ts.day.takingoff(g.day>datetime('15-July-2018') & ~isnan(Ts.day.takingoff'));
disp(['50% of all bird taking-off in Autumn happened in ' num2str(sum( (cumsum(sort(tmp,'descend'),'omitnan') ./ nansum(tmp)) < .5) )])
Accumulation of bird in Spring
spring_entering = nansum(Ts.Fout.day.entering(g.day<datetime('15-July-2018'))/1000000)
spring_leaving = nansum(Ts.Fout.day.leaving(g.day<datetime('15-July-2018'))/1000000)
spring_delta = spring_entering + spring_leaving
Accumulation (or dissipation) of bird in Autumns
autumn_entering = nansum(Ts.Fout.day.entering(g.day>datetime('15-July-2018'))/1000000)
autumn_leaving = nansum(Ts.Fout.day.leaving(g.day>datetime('15-July-2018'))/1000000)
autumn_delta = autumn_entering + autumn_leaving
Recrutement = Flux out/Flux in = Bird in + reproduction / Bird in
disp(['Ratio of Aumtumn departure over spring arrival: ' num2str(mean(-autumn_leaving ./ spring_entering)) ])
disp(['Ratio of Aumtumn accumulation over spring accumulation: ' num2str(mean(-autumn_delta ./ spring_delta))])
Seasonal Flow: In and Out
Flux per sector (15min)
transect_name={'North', 'East', 'Alpes', 'Spain' ,'Altlantic','England'};
transect(:,:,1)= [42 51 ; 43 82];
transect(:,:,2)= [20 45 ; 78 87];
transect(:,:,3)= [1 19 ; 40 84];
transect(:,:,4)= [1 9; 1 39];
transect(:,:,5) = [5 24; 1 18];
transect(:,:,6)= [25 43 ; 1 45];
sec_col={'r','b','y','g','c','k'};
cat=zeros(gext.nlat,gext.nlon);
for i_s=1:numel(transect_name)
cat(transect(1,1,i_s):transect(1,2,i_s),transect(2,1,i_s):transect(2,2,i_s))=i_s;
end
Fout.sec=nan(g.nat,numel(transect_name),2);
for i_s=1:numel(transect_name)
Fout.sec(:,i_s,1) = sum(reshape(Fout.day.entering(repmat(cat_trans==i_s,1,1,g.nat)),[],g.nat));
Fout.sec(:,i_s,2) = sum(reshape(Fout.day.leaving(repmat(cat_trans==i_s,1,1,g.nat)),[],g.nat));
end
Fout.sec(~knownday,:,:)=nan;
Flux per sector (season)
Fout.sec_season=nan(numel(transect_name),4);
Fout.sec_season(:,1)=nansum(Fout.sec(g.day<datetime('2018-07-01'),:,1));
Fout.sec_season(:,2)=nansum(Fout.sec(g.day<datetime('2018-07-01'),:,2));
Fout.sec_season(:,3)=nansum(Fout.sec(g.day>datetime('2018-07-01'),:,1));
Fout.sec_season(:,4)=nansum(Fout.sec(g.day>datetime('2018-07-01'),:,2));
Figure illustration sector
figure('position',[0 0 1000 400]);
h = worldmap([g.lat(1)-10 g.lat(end)+10], [g.lon(1)-10 g.lon(end)+10]);
%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')
figure; worldmap([g.lat(1)-1 g.lat(end)+1], [g.lon(1)-1 g.lon(end)+2]);
tmp = cat; tmp(mean(Fout.Fout>0,3)<.02)=nan;
surfm([gext.lat2D gext.lat; 55.5*ones(1,gext.nlon+1)]-0.125, [gext.lon2D 14.5*ones(gext.nlat,1); gext.lon' 14.5]-0.125 ,padarray(tmp,[1 1],nan,'post'),'Facecolor','flat');
c=colorbar; c.Ticks=1:6; c.TickLabels=transect_name;
Flux per sector (daily)
figure('position',[0 0 1000 1000]);
for i_s=1:numel(transect_name)
subplot(numel(transect_name),1,i_s); hold on;
bar(g.day, Fout.sec(:,i_s,1)/1000000,'FaceColor',col2(8,:));
bar(g.day, -Fout.sec(:,i_s,2)/1000000,'FaceColor',col2(10,:));
ylabel(transect_name{i_s})
ylim([0 30])
box on; grid on
end