A simple method to separate birds and insects in single-pol weather radar

Raphaël Nussbaumer1,*, Baptiste Schmid1 , Silke Bauer1 and Felix Liechti1

1 Swiss Ornithological Institute, Sempach, Switzerland

* Correspondence: raphael.nussbaumer@vogelwarte.ch

Introduction

This code accompagnes the publication: Nussbaumer, R.; Schmid, B.; Bauer, S.; Liechti, F. Technical a Gaussian Mixture Model to Separate Birds and Insects in Single-Polarization Weather Radar Data. Remote Sens. 2021, 13, 1989. https://doi.org/10.3390/rs13101989

It explains how the methodology is implemented, how the figure are produced and how the resulting dateset is exported (available on zenodo: https://doi.org/10.5281/zenodo.3610184)

Abstract

Recent and archived data from weather radar networks are extensively used for the quantification of continent-wide bird migration patterns. While the process of discriminating birds from weather signals is well established, insect contamination is still a problem. We present a simple method combining two Doppler radar products within a Gaussian mixture model to estimate the proportions of birds and insects within a single measurement volume, as well as the density and speed of birds and insects. This method can be applied to any existing archives of vertical bird profiles, such as the European Network for the Radar surveillance of Animal Movement repository, with no need to recalculate the huge amount of original polar volume data, which often are not available.

Table of Contents

Introduction Abstract Setup Load the data Determine windspeed at radar location Compute airspeed and standar deviation v_a radial velocity sdvvp Model fitting Compute empirical PDF and fit with a mixture of Gaussian Compute amplitude ratio Compute amplitude ratio over time Compute amplitude ratio over space Compute amplitude ratio over altitude Compute amplitude ratio per radar and month Temporal interpolation of the amplitude ratio Separation of insect and weather Correction for bird and insect speed Apply model to the dataset Export the data Result

Setup

Load the data

load('data/dc_corr')
load('coastlines.mat')

Define default color as viridis. Code can be downloadeded from https://www.mathworks.com/matlabcentral/fileexchange/51986-perceptually-uniform-colormaps

if (exist('viridis'))
    set(0,'DefaultFigureColormap',viridis);
else
    web('https://www.mathworks.com/matlabcentral/fileexchange/51986-perceptually-uniform-colormaps')
end

Determine windspeed at radar location

The U (parmID=131) and V (parmID=132) component of wind are downloaed from the ERA5 reanalysis (https://doi.org/10.24381/cds.bd0915c6) at pressure level (from 1000hPa to 550hPa), downloaded at the maximal resoluation (hourly and 0.25°x0.25°). See python file ECMWF/script_ECMWF.py for more details on the download.

%ncdisp(file{1});
file={'./ECMWF/2018_pressure_1.nc','./ECMWF/2018_pressure_2.nc','./ECMWF/2018_pressure_3.nc','./ECMWF/2018_pressure_4.nc'};
wind.time = datenum('01-janv-2018'):1/24:datenum('31-dec-2018 23:00');
wind.latitude=flip(double(ncread(file{1},'latitude')));
wind.longitude=double(ncread(file{1},'longitude'));
wind.pressure=double([ncread(file{1},'level') ; ncread(file{2},'level') ; ncread(file{3},'level') ; ncread(file{4},'level') ]);
wind.alt = (1-(wind.pressure*100/101325).^(1/5.25588))/2.25577/10^(-5);
tmp_1 = permute(flip(ncread(file{1},'u'),2) , [2 1 4 3]);
tmp_2 = permute(flip(ncread(file{2},'u'),2) , [2 1 4 3]);
tmp_3 = permute(flip(ncread(file{3},'u'),2) , [2 1 4 3]);
tmp_4 = permute(flip(ncread(file{4},'u'),2) , [2 1 4 3]);
wind.u = cat(4,tmp_1,tmp_2,tmp_3,tmp_4); % m/s
tmp_1 = permute(flip(ncread(file{1},'v'),2) , [2 1 4 3]);
tmp_2 = permute(flip(ncread(file{2},'v'),2) , [2 1 4 3]);
tmp_3 = permute(flip(ncread(file{3},'v'),2) , [2 1 4 3]);
tmp_4 = permute(flip(ncread(file{4},'v'),2) , [2 1 4 3]);
wind.v = cat(4,tmp_1,tmp_2,tmp_3,tmp_4); % m/s
tmp_1 = permute(flip(ncread(file{1},'t'),2) , [2 1 4 3]);
tmp_2 = permute(flip(ncread(file{2},'t'),2) , [2 1 4 3]);
tmp_3 = permute(flip(ncread(file{3},'t'),2) , [2 1 4 3]);
tmp_4 = permute(flip(ncread(file{4},'t'),2) , [2 1 4 3]);
wind.t = cat(4,tmp_1,tmp_2,tmp_3,tmp_4); % K
clear tmp_* file

All three variables are linearly interpolated (time-space 4D) at each datapoint of the weather radar data.

Fu = griddedInterpolant({wind.latitude,wind.longitude,datenum(wind.time),wind.alt},wind.u,'linear','linear');
Fv = griddedInterpolant({wind.latitude,wind.longitude,datenum(wind.time),wind.alt},wind.v,'linear','linear');
Ft = griddedInterpolant({wind.latitude,wind.longitude,datenum(wind.time),wind.alt},wind.t,'linear','linear');
for i_d=1:numel(dc)
    windu = permute(Fu({dc(i_d).lat,dc(i_d).lon,datenum(dc(i_d).time),dc(1).alt}),[3,4,1,2]);
    windv = permute(Fv({dc(i_d).lat,dc(i_d).lon,datenum(dc(i_d).time),dc(1).alt}),[3,4,1,2]);
    dc(i_d).wt = permute(Ft({dc(i_d).lat,dc(i_d).lon,datenum(dc(i_d).time),dc(1).alt}),[3,4,1,2]);
    dc(i_d).ws = windu + 1i*windv;
end

Export dataset and clear

clear wind Fu Fv
% save('data/dc_corr','dc','start_date','end_date','quantity','-v7.3')

Compute airspeed and standar deviation v_a radial velocity sdvvp

v_a = nan(numel(dc(1).time), numel(dc(1).alt), numel(dc));
sdvvp = v_a;
for i_d=1:numel(dc)
    
    % Compute the n/s and e/w componenent of the flight speed
    % [-u,v] = pol2cart(deg2rad(dc(i_d).dd+90), dc(i_d).ff);
    % u =  dc(i_d).ff .* sind(dc(i_d).dd); % m/s | 0° is north and 90° is east. -> u is east (+) - west (-)
    % v = dc(i_d).ff .* cosd(dc(i_d).dd); % m/s  -> v is north (+) - south (-)
    % gs = u + 1i*v;
    gs =  dc(i_d).ff .* exp(deg2rad(mod(-dc(i_d).dd+90,360))*1i);
    
    % Remove data not cleaned for density (e.g. rain).
    id = isnan(dc(i_d).dens3) | dc(i_d).dens3==0;
    gs(id) = nan;
    sd_vvp = dc(i_d).sd_vvp;
    sd_vvp(id) = nan;
    
    v_a(:,:,i_d) = gs - dc(i_d).ws;
    sdvvp(:,:,i_d) = sd_vvp;
end

Plot airpseed and sdvvp histogram

figure('position',[0 0 800 400]); hold on;

histogram(abs(v_a),'EdgeColor','none','DisplayName','Airspeed')

histogram(sdvvp,'EdgeColor','none','DisplayName','SD radial velocity')

legend; xlabel('speed [m/s]'); xlim([0 25])

Model fitting

Compute empirical PDF and fit with a mixture of Gaussian

Avoid to recompute the gmfit and the ampli* and f*

load('data/insect_removal')

Define variable used later

t =  datenum(dc(1).time)-datenum(dc(1).time(1));
nmonth = 12;
tmonth = linspace(t(1), t(end), nmonth+1);
tweek = t(1):7:t(end);
nweek = numel(tweek);
tmonth_mid = tmonth(1:end-1)+diff(tmonth)/2;
tmonth_mid_rep = [tmonth_mid-365 tmonth_mid tmonth_mid+365]

France and German radar

id_de = [false; true(15,1); false(21,1)];
id_fr = [false(16,1); true(19,1); false(2,1)];

Build the vector of data removing nan value

X = [abs(v_a(:)) sdvvp(:)];
X = X(~any(isnan(X),2),:);

Fit a kernel density function

[xi1, xi2] = meshgrid(0:.25:15,0:.1:7);
ft = reshape(ksdensity(X, [xi1(:), xi2(:)]),size(xi1));
ft = ft./sum(ft(:));

Fit Gaussian mixture model with two components

S.mu = [8 4 ; 2.5 3];
S.Sigma = cat(3,[11 -1 ; -1 1], [2 0.1 ; 0.1 1]);
S.ComponentProportion = [0.7 0.3];
rng('default')
gmfit = fitgmdist(X,2,'Replicates',1,'Start',S);

Compute the normalized pdf of each gaussian

gm1 = gmdistribution(gmfit.mu(1,:), gmfit.Sigma(:,:,1));
f_1 = reshape(pdf(gm1,[xi1(:) xi2(:)]),size(xi1));
f_1 = f_1./sum(f_1(:));
gm2 = gmdistribution(gmfit.mu(2,:), gmfit.Sigma(:,:,2));
f_2 = reshape(pdf(gm2,[xi1(:) xi2(:)]),size(xi1));
f_2 = f_2./sum(f_2(:));

Plot the fitted pdf

figure('position',[0 0 300 300]); hold on; hold on

% plot(xi1(islocalmax(ft)&islocalmax(ft,2)), xi2(islocalmax(ft)&islocalmax(ft,2)),'.r')

surf(xi1, xi2, ft,'FaceAlpha',0.5)

[~,p1]=contour3(xi1, xi2, gmfit.ComponentProportion(1)*f_1,10,'Color',[162 29 49]/255,'linewidth',1.5);

[~,id_max] = max(f_1(:));

plot3(xi1(id_max),xi2(id_max),gmfit.ComponentProportion(1)*f_1(id_max),'+','Color',[162 29 49]/255,'LineWidth',2,'MarkerSize',6)

plot3([xi1(id_max) xi1(id_max)],[xi2(id_max) xi2(id_max)], [0 gmfit.ComponentProportion(1)*f_1(id_max)],'Color',[162 29 49]/255,'LineWidth',1.5 )

[~,p2]=contour3(xi1, xi2, gmfit.ComponentProportion(2)*f_2,10,'Color',[127 47 141]/255,'linewidth',1.5);

[~,id_max] = max(f_2(:));

plot3(xi1(id_max),xi2(id_max),gmfit.ComponentProportion(2)*f_2(id_max),'+','Color',[127 47 141]/255,'LineWidth',2,'MarkerSize',6)

plot3([xi1(id_max) xi1(id_max)],[xi2(id_max) xi2(id_max)], [0 gmfit.ComponentProportion(2)*f_2(id_max)],'Color',[127 47 141]/255,'LineWidth',1.5 )

axis tight;shading interp; view(37,31)

xlabel('airspeed [m/s]'); ylabel('\sigma_{vvp} [m/s]'); zlabel('PDF')

%legend([p1 p2], 'Multi-normal corresponding to bird','Multi-normal corresponding to insect')

Ax = gca; Ax.ZAxis.TickValues=[];

colormap(gca,viridis) %colormap(gca,1-(1-viridis)/1.2) %

Display parameters of the gaussian

disp(gmfit.mu)
    7.9599    4.1032
    2.6119    2.7799
disp(gmfit.Sigma)
(:,:,1) =

   11.5652   -1.1574
   -1.1574    0.9159


(:,:,2) =

    1.8133    0.1578
    0.1578    1.0829

Plot figure as a comparison with threashold

figure('position',[0 0 800 600]); hold on; hold on

h(1)=pcolor(xi1(1,:), xi2(:,1), ft./max(ft(:)));

h(2)=plot([0 15],[2 2],'r','linewidth',2);

plot([5 5],[0 7],'r','linewidth',2);

[~,h(3)]=contour(xi1, xi2, f_1,6,'Color',[162 29 49]/255,'LineWidth',1.5);

[~,id_max] = max(f_1(:)); plot(xi1(id_max),xi2(id_max),'+','Color',[162 29 49]/255,'LineWidth',2,'MarkerSize',10)

[~,h(4)]=contour(xi1, xi2, f_2,6,'Color',[127 47 141]/255,'LineWidth',1.5);

[~,id_max] = max(f_2(:)); plot(xi1(id_max),xi2(id_max),'+','Color',[127 47 141]/255,'LineWidth',2,'MarkerSize',10)

[C,h(5)]=contour(xi1, xi2, f_2./(f_1+f_2),.1:.2:.9,'w','ShowText','on','LineWidth',2);

clabel(C,h(5),'FontSize',15,'Color','w','LabelSpacing',500,'FontWeight','bold')

axis tight;shading interp;

xlabel('airspeed [m/s]'); ylabel('\sigma_{vvp} [m/s]');

legend(h,{'Kernel density estimation','Traditional Threashold approach','Bird Gaussian Fit','Insect Gaussian Fit','New Proportion approach'})

colormap(gca,1-(1-viridis)/1.2) %colormap(gca,viridis)

Compute amplitude ratio

Compute amplitude ratio over time

amplit = nan(1,nmonth);
ftime = nan(size(xi1,1), size(xi1,2),nmonth);
for i=1:nmonth
    id = tmonth(i)<t & tmonth(i+1)>t;
    X = [reshape(abs(v_a(id,:,:)),[],1),reshape(sdvvp(id,:,:),[],1)];
    X = X(~any(isnan(X),2),:);
    if ~isempty(X)
        ftime(:,:,i) = reshape(ksdensity(X, [xi1(:), xi2(:)]),size(xi1));
        ftime(:,:,i) = ftime(:,:,i)./sum(sum(ftime(:,:,i)));
        mismatch = @(alpha) ( alpha.*f_1 + (1-alpha).*f_2);
        amplit(i) = fminsearch(@(alpha) sum(sum((ftime(:,:,i)-mismatch(alpha)).^2)) , 0.5 );
    end
end
amplit(amplit==0)=nan;

figure('position',[0 0 1200 800]); hold on;
c_map = magma;
for i=1:nmonth
    ax=subplot(3,4,i); hold on
    if i==1
        colorbar
        colormap(gca,magma)
    else
        id = tmonth(i)<t & tmonth(i+1)>t;
        pcolor(xi1(1,:), xi2(:,1), ftime(:,:,i)./max(max(ftime(:,:,i))))
        contour(xi1, xi2, amplit(i)*f_1,[0.0002 0.01],'Color',[162 29 49]/255,'LineWidth',1.5);
        [~,id_max] = max(f_1(:));
        plot(xi1(id_max),xi2(id_max),'+','Color',[162 29 49]/255,'LineWidth',2,'MarkerSize',amplit(i)*20)
        contour(xi1, xi2, (1-amplit(i))*f_2,[0.0002 0.01],'Color',[127 47 141]/255,'LineWidth',1.5);
        [~,id_max] = max(f_2(:));
        plot(xi1(id_max),xi2(id_max),'+','Color',[127 47 141]/255,'LineWidth',2,'MarkerSize',(1-amplit(i))*20)
        [C,h]=contour(xi1, xi2, (1-amplit(i))*f_2./(amplit(i)*f_1+(1-amplit(i))*f_2),[.1 .5 .9],'w','ShowText','on','LineWidth',2);
        clabel(C,h,'Color','w','LabelSpacing',500,'FontWeight','bold')
        axis tight;shading interp;
        title(month(median(dc(1).time(id)),'name'))
        %xlabel('airspeed [m/s]'); ylabel('\sigma_{vvp} [m/s]');
        colormap(gca,1-(1-viridis)/1.2)
        
        box on
        i_c = ceil(amplit(i)*size(c_map,1));
        ax.XColor = c_map(i_c,:);
        ax.YColor = c_map(i_c,:);
        ax.LineWidth=10;
        ax.XTick=[]; ax.YTick=[];
    end
end