Inference of Model

Load and import dataset

clear all; load('./data/dc_corr.mat','dc'); addpath('./functions/'); load('coastlines.mat'); warning('off')
load('./data/BelowRadarMPS.mat');
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) );

Convert the data to a simple structure

Define the data structure with all data.
data.name = {dc.name}';
data.lat = [dc.lat]';
data.lon = [dc.lon]';
data.height = [dc.height]';
data.heightDEM = [dc.heightDEM]';
data.maxrange = [dc.maxrange]';
data.nrad = numel(data.name);
data.time = dc(1).time';
data.alt = dc(1).alt';
data.NNT = datenum(repmat(data.time,1,numel(dc))-mean(cat(3,[dc.dawn],[dc.dusk]),3)) ./ datenum([dc.dawn]-[dc.dusk])*2;
Initiate empty density matrix
data_denss_sim = nan(numel(data.time), numel(dc), size(MPS{1},3));
Loop for eac radar
for i_d=1:numel(dc)
% Vertical agregation accounts for the DEM based on VolBelow.
% First, find the direction of flight of the lowest non-nan bin
[~, id] = max(~isnan(dc(i_d).dd), [], 2);
dir_lowest = dc(i_d).dd(sub2ind(size(dc(i_d).dd), (1:numel(data.time))',id));
% Transform this direction angle (0-360°) into the index of alpha (-90°-90°)
dir_lowest_round = round(dir_lowest,-1);
id_alpha=nan(size(dir_lowest_round));
id_alpha(dir_lowest_round>=0 & dir_lowest_round<90) = dir_lowest_round(dir_lowest_round>=0 & dir_lowest_round<90)+90;
id_alpha(dir_lowest_round>=90 & dir_lowest_round<270) = dir_lowest_round(dir_lowest_round>=90 & dir_lowest_round<270)-90;
id_alpha(dir_lowest_round>=270 & dir_lowest_round<360) = dir_lowest_round(dir_lowest_round>=270 & dir_lowest_round<360)-270;
id_alpha = id_alpha/10+1;
height_vol = nan(numel(data.time),numel(data.alt));
height_vol(~isnan(id_alpha),:) = dc(i_d).VolBelow(id_alpha(~isnan(id_alpha)),:);
% Integrate volume density into surface density based on the direction and corresponding volume. bird/km^2
% data.denss(:,i_d) = sum(data_dens(:,:,i_d) .* height_vol,2);
% Computing with result of MPS
tmp = dc(i_d).dens4;
data_denss_MPS = nan(size(tmp,1),size(MPS{i_d},3),size(tmp,2));
for i_real=1:size(MPS{i_d},3)
tmp(:,1:dc(i_d).scatter_lim-1) = MPS{i_d}(:,1:dc(i_d).scatter_lim-1,i_real);
data_denss_MPS(:,i_real,:) = tmp .* height_vol;
end
data_denss_sim(:,i_d,:) = sum(data_denss_MPS,3);
% Compute speed
nb_bird = permute(nanmean(data_denss_MPS,2),[1 3 2]);
w_MTR = nb_bird ./ repmat(nansum(nb_bird,2),1,size(nb_bird,2));
data.vs(:,i_d) = nansum(w_MTR .* dc(i_d).v2,2);
data.us(:,i_d) = nansum(w_MTR .* dc(i_d).u2,2);
% if data present initially
data.isdata(:,i_d)=~all(isnan(dc(i_d).DBZH),2);
end
data.vs(data.vs==0)=nan;
data.us(data.us==0)=nan;
Set data during the day to zero
data_denss_sim(repmat((data.NNT<-1 | data.NNT>1),1,1,size(data_denss_sim,3)))=nan;
data.vs(data.NNT<-1 | data.NNT>1)=nan;
data.us(data.NNT<-1 | data.NNT>1)=nan;
Remove 0 value because of the transformation
data_denss_sim(data_denss_sim==0)=nan;
Remove data when there is only a few data point (8 => 2h00) per day
tmp=dateshift(data.time,'start','day','nearest');
day_id=findgroups(datenum(tmp));
G = findgroups(day_id);
nb = splitapply(@nansum,~isnan(data_denss_sim),G);
data_denss_sim(nb(day_id,:,:)<8)=nan;
Remove time with all day or no data to save space
data.id_t=all(all(isnan(data_denss_sim),2),3);
data.time(data.id_t)=[];
data.ntime = numel(data.time);
data.NNT(data.id_t,:)=[];
data_denss_sim(data.id_t,:,:)=[];
data.vs(data.id_t,:)=[];
data.us(data.id_t,:)=[];
data.isdata(data.id_t,:)=[];
tmp=dateshift(data.time,'start','day','nearest');
[data.day_id,data.day]=findgroups(datenum(tmp));
G = findgroups(data.day_id);
nb = splitapply(@nansum,~isnan(data_denss_sim),G);
Clear variable to make space (if necessary)
clear id_t tmp1 dens_lim tmp w_MTR nb_bird data_u data_v ff dd data_dens dir_lowest dir_lowest_round height_vol id_alpha id i_d dir i_real data_denss_MPS MPS

Rain

Determine the threashold
file='./ECMWF/2018_srf.nc'; %ncdisp(file);
rain.time = datenum('01-janv-2018'):1/24:datenum('31-dec-2018 23:00');
rain.latitude=flip(double(ncread(file,'latitude')));
rain.longitude=double(ncread(file,'longitude'));
rain.tp = permute(flip(ncread(file,'tp'),2) , [2 1 3]); % Total precipitation [m] (accumulation over the hour ending at the forcast step)
F = griddedInterpolant({rain.latitude,rain.longitude,datenum(rain.time)},rain.tp);
for i_d=1:numel(dc)
data.rain(:,i_d) = F({data.lat(i_d),data.lon(i_d),datenum(data.time)});
end
Consider only night time, where data is present (i.e. DBZH available)
consider = data.NNT>-1 & data.NNT<1 & data.isdata;
Bird density
tmp = data_denss_sim(:,:,1);
Figure to illustrate
% figure; hold on;
% histogram(data.rain(isnan(tmp)&consider))
% histogram(data.rain(~isnan(tmp)&consider))
% set(gca,'YScale','log')
% figure; scatter(data.rain(consider), tmp(consider),'k.')
figure('position',[0 0 1000 400]); hold on;
plot(tmp(:,2))
plot(data.isdata(:,2)*100)
yyaxis right;
plot(data.rain(:,2))
thr=0:0.00001:0.005;
score=nan(3,numel(thr));
for i=1:numel(thr)
rain_id=data.rain>=thr(i);
score(1,i)=mean(isnan(tmp(rain_id&consider)));
score(2,i)=mean(isnan(tmp(~rain_id&consider)));
score(3,i)=sum(rain_id(:)&consider(:))/sum(consider(:));
end
figure('position',[0 0 1000 400]); plot(thr,score'); xlabel('threshold value [m/hr]'); ylabel('%')
legend('% of data nan amongst the eliminated data','% of nan in the non-eliminated data')
This threashold value is defined at 1 mm/hr -> 0.001 m at the hourly resolution of the total precipitation according to https://confluence.ecmwf.int//display/CKB/ERA5+data+documentation#ERA5datadocumentation-Meanratesandaccumulations
data.mask_rain_thr=0.001;
Clear variable
clear rain F file consider tmp score

Transformation of data

Raw data very skew, transformation to normal space using a kernel-based density fit.
Transformation does not work for dens=0. remove for the fit
pd=fitdist(data_denss_sim(:),'Kernel','Support','positive');
Too many data for the forward transformation. Reduce the site of the data to transform using only the value rounded to 2 decimal. Compute the transformation of the MPS simulation and compute the std
trans.dens_axis = logspace(log10(min(data_denss_sim(:))), log10(max(data_denss_sim(:))),1000);
trans.denst_axis = norminv(cdf(pd,trans.dens_axis));
trans.dens_axis([false diff(trans.denst_axis)==0])=[];
trans.denst_axis([false diff(trans.denst_axis)==0])=[];
data.denst = interp1(trans.dens_axis,trans.denst_axis,data_denss_sim,'pchip');
data.denst_m = nanmean(data.denst,3);
% data.denst_s = nanstd(data.denst,[],3);
inverse: SUPER SLOW
% data.dens_m = icdf(pd,normcdf(data.denst_m));
data.dens_m = interp1(trans.denst_axis,trans.dens_axis,data.denst_m,'pchip');
Figure
figure('position',[0 0 1000 400]); hold on;
subplot(1,2,1); hold on;
histogram(data_denss_sim(:),'Normalization','probability'); plot(0:500,pdf(pd,0:500))
tmp = nanmean(data.denst,3);
subplot(1,2,2); histfit(tmp(:),100)
clear pd

Block decompostion

We compute the parameters of the trends and covariance function on each temporal subset (=block) of the dataset. Each block of the dataset is of size data_block days and add some overlap with previous and next day of block_overlap days.
data.block.date = [737121 737184 737273 737384]; % '13-Feb-2018', '02-Mar-2018', '23-Apr-2018','24-Sep-2018', '20-Nov-2018', '01-Jan-2019'
data.block.date_str={'Spring (02-Mar - 04-May)','Summer (4-May - 1-Aug)','Autumn (1-Aug - 20-Nov)','Winter (20-Nov - 02-Mar)'};
data.block.col = [[81 161 114];[206 202 105];[213 122 93];[85 126 170]]/255;
These date are defined with the figure below
U=nanmean(splitapply(@nanmean, (data.us*60*60/1000.*data.dens_m), data.day_id),2);
V=nanmean(splitapply(@nanmean, (data.vs*60*60/1000.*data.dens_m), data.day_id),2);
MTR = sqrt(U.^2+V.^2);
figure('position',[0 0 1000 400]); hold on
% Background color
for i_b=1:numel(data.block.date)-1
fill([data.block.date(i_b)-.5 data.block.date(i_b)-.5 data.block.date(i_b+1)-.5 data.block.date(i_b+1)-.5],[4.5 0 0 4.5]*1000,data.block.col(i_b,:))
end
fill([min(data.day)-.5 min(data.day)-.5 data.block.date(1) data.block.date(1)-.5],[4.5 0 0 4.5]*1000,data.block.col(4,:))
fill([data.block.date(end)-.5 data.block.date(end)-.5 max(data.day) max(data.day)-.5],[4.5 0 0 4.5]*1000,data.block.col(4,:))
% MTR
bar(data.day,MTR,'k')
ylabel('Average movement of bird per day [bird/km/hr] (MTR)');
% MTR direction
yyaxis right; plot(data.day,atan2d(V,U),'.','Color',[.2 .2 .2]);
ylabel('Direction of the movement (0°)');
line([data.day(1) data.day(end)],[0 0],'Color',[.2 .2 .2]); text(data.day(1),8,'East','Color',[.2 .2 .2])
line([data.day(1) data.day(end)],[90 90],'Color',[.2 .2 .2]); text(data.day(1),98,'North','Color',[.2 .2 .2])
line([data.day(1) data.day(end)],[-180 -180],'Color',[.2 .2 .2]); text(data.day(1),-172,'West','Color',[.2 .2 .2])
line([data.day(1) data.day(end)],[-90 -90],'Color',[.2 .2 .2]); text(data.day(1),-82,'South','Color',[.2 .2 .2])
ylim([-180 180]);
datetick('x'); xlabel('Date'); ax = gca;ax.YAxis(1).Color = 'k';ax.YAxis(2).Color = 'k'; axis tight;
for i_b=1:numel(data.block.date)
if i_b==numel(data.block.date)
data.b(datenum(data.time)>=data.block.date(4) | datenum(data.time)<data.block.date(1)) = i_b;
else
data.b(datenum(data.time)>=data.block.date(i_b) & datenum(data.time)<data.block.date(i_b+1)) = i_b;
end
end

Average Daily

Compute the radar and day index (t and r)
multi.t = data.day;
multi.r = 1:data.nrad;
[multi.R,multi.T]=meshgrid(multi.r,multi.t);
Compute the daily average
G = findgroups(data.day_id);
multi.M = splitapply(@nanmean,data.denst,G);
nb = splitapply(@nansum,~isnan(data.denst),G);
multi.M_m = mean(multi.M,3);
multi.M_var = var(multi.M,[],3);
multi.isnan = isnan(multi.M_m);
Figure
figure('position',[0 0 1000 400]); hold on;
subplot(2,1,1); imagesc(multi.M_m')
subplot(2,1,2); imagesc(sqrt(multi.M_var)'); caxis([0 .1])
Attribute block
multi.b = nan(size(multi.t));
for i_b=1:numel(data.block.date)
if i_b==numel(data.block.date)
multi.b(multi.t>=data.block.date(i_b) | multi.t<data.block.date(1)) = i_b;
else
multi.b(multi.t>=data.block.date(i_b) & multi.t<data.block.date(i_b+1)) = i_b;
end
end
multi.B = repmat(multi.b,1,data.nrad);
data.day_b = multi.b;

Trend

The trend is commputed as a polynomila of degree poly_t_degree in the time axis and poly_t_degree in space (2D).
multi_poly_t_degree=0;
multi_poly_d_degree=1;
The trend is determined with a Generalized Linear Square regression (GLS), where the following linear problem is solved
,
where are the coefficient of the polynomes, C is the covariance matrix of the data y. Xis the Vandermonde matrix for the time and coordinate of the data y. In our case:
X = trend(multi.T,data.lat(multi.R),data.lon(multi.R),multi_poly_t_degree,multi_poly_d_degree);
C = diag(1+multi.M_var(:));
y = multi.M_m(:);
We can then create the forward function of the trend which given any location (time and coordinate) and polynomial coefficinet w return the value of the trend
multi.f_trend = @(t,lat,lon,beta) reshape(trend(t(:),lat(:),lon(:),multi_poly_t_degree,multi_poly_d_degree) * beta(:), size(t));
We will then compute the detrended value of the multi-night scale componenent
multi.M_mn=nan(size(multi.M_m));
multi.beta = nan(multi_poly_t_degree+ (multi_poly_d_degree)*2+1,numel(data.block.date));
idm_b=cell(numel(data.block.date),1);
We loop over all block...
for i_b=1:numel(data.block.date)
... first compute the block index....
id = find(~multi.isnan & multi.B==i_b);
.... and then, the trend coefficient
multi.beta(:,i_b) = ((X(id,:)'*(C(id,id)\X(id,:)))\X(id,:)')*(C(id,id)\y(id)); % GLS
Figure
if false
figure('position',[0 0 1000 400]); hold on;
t_tmp = sort(unique(multi.T(id)));
plot(t_tmp,multi.f_trend(t_tmp,mean(data.lat)*ones(numel(t_tmp),1),mean(data.lon)*ones(numel(t_tmp),1), multi.beta(:,i_b)),'-r')
plot(multi.T(id),y(id),'.g')
plot(multi.T(id),y(id),'.k')
figure('position',[0 0 1000 400]);
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')
[lat_tmp,lon_tmp]=meshgrid(floor(min(data.lat)):.2:ceil(max(data.lat)), floor(min(data.lon)):.2:ceil(max(data.lon)));
S=multi.f_trend(mean(multi.T(id))*ones(size(lat_tmp)),lat_tmp,lon_tmp,multi.beta(:,i_b));
surfm(lat_tmp,lon_tmp,S)
plotm(coastlat, coastlon,'k'); geoshow('worldrivers.shp','Color', 'blue')
[G,id]=findgroups(multi.R(id));
scatterm(data.lat(id),data.lon(id),100,splitapply(@mean,y(id), G),'filled','MarkerEdgeColor','k');
caxis([min(y(id)) max(y(id))])
end
Compute the de-trended value
multi.M_mn(id) = multi.M_m(id)-multi.f_trend(multi.T(id),data.lat(multi.R(id)),data.lon(multi.R(id)),multi.beta(:,i_b));
end
Figure
% col = [[0.0235 0.4824 0.2157];[0.7333 0.7059 0.1569];[0.7725 0.2549 0.1020];[0.0510 0.2784 0.5255]];
figure('position',[0 0 1000 400]); 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)],[4 -4 -4 4],data.block.col(i_b,:))
end
fill([min(data.day) min(data.day) data.block.date(1) data.block.date(1)],[4 -4 -4 4],data.block.col(4,:))
fill([data.block.date(end) data.block.date(end) max(data.day) max(data.day) ],[4 -4 -4 4],data.block.col(4,:))
plot(multi.T(:), y,'.k')
for i_b=1:numel(data.block.date)
id = find(~multi.isnan & multi.B==i_b);
t_tmp = sort(unique(multi.T(id)));
plot(t_tmp, multi.f_trend(t_tmp,mean(data.lat(multi.R(id)))*ones(numel(t_tmp),1),mean(data.lon(multi.R(id)))*ones(numel(t_tmp),1), multi.beta(:,i_b)),'-r','Linewidth',2)
end
datetick('x'); ylim([-4 4]); 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(~multi.isnan & multi.B==i_b);
subplot(2,2,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,multi.f_trend(mean(multi.T(id))*ones(size(LAT)),LAT,LON,multi.beta(:,i_b)))
plotm(coastlat, coastlon,'k'); geoshow('worldrivers.shp','Color', 'blue')
[G,tmp]=findgroups(multi.R(id));
tmp2 = splitapply(@mean,multi.M(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]); imagesc(multi.M_mn'); xlabel('Time'); ylabel('Radar')