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

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