Birds at risk

Determine how many birds are at risk of collision with windturbine where and when.
Table of Contents

Setup

Load data

Load windfarm data which has been processed (windfram_processing.mlx)
load('data/windfarms_processed.mat')
load('data/energy_processed.mat');
Load the grid of the bird density estimation
load('../2018/data/Density_estimationMap','g');
Load flow estimation (density and speed)
load('../2018/data/Density_estimationMap','gd');
load('../2018/data/Flight_estimationMap','guv');
Load the ratio map of altiudinal distribution of the bird
load('./data/ratio');
Add path
addpath('functions')

Build other usefull data for later

define colormap
greens = interp1([1 128 256],[220, 233, 187;28, 155, 117;2, 76, 53]/255,1:256);
reds = interp1([1 128 256],[239, 220, 217;181, 57, 104;119, 38, 76]/255,1:256);
blues = interp1([1 128 256],[197, 221, 244;18, 122, 161;45, 66, 88]/255,1:256);
Compute area
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);

Create country map

bd = load('borderdata.mat');
inCountry=false(numel(g.lat),numel(g.lon),3);
ind = strcmp(bd.places,'France');
inCountry(:,:,1) = reshape(inpolygon(g.lon2D(:), g.lat2D(:), bd.lon{ind}', bd.lat{ind}'),size(g.lon2D));
ind = strcmp(bd.places,'Germany');
inCountry(:,:,2) = reshape(inpolygon(g.lon2D(:), g.lat2D(:), bd.lon{ind}', bd.lat{ind}'),size(g.lon2D));
ind = strcmp(bd.places,{'Belgium'}) | strcmp(bd.places,{'Netherlands'}) | strcmp(bd.places,{'Luxembourg'});
inCountry(:,:,3) = reshape(inpolygon(g.lon2D(:), g.lat2D(:), [bd.lon{ind}]', [bd.lat{ind}]'),size(g.lon2D));
Illustration
figure('position',[0 0 500 500]); hold on;
h = worldmap([g.lat(1)-2 g.lat(end)+2], [g.lon(1)-2 g.lon(end)+2]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
bordersm('countries','facecolor',[228 238 243]./255);
surfm(g.lat,g.lon,sum(inCountry.*reshape(1:3,1,1,[]),3))
bordersm('countries','k');

Compute distance to coastline

load coastlines
dcoast = pdist2([g.lat2D(:) g.lon2D(:)],[coastlat coastlon], @lldistkm);
dcoast = reshape(min(dcoast,[],2),size(g.lat2D));
Illustration
figure('position',[0 0 500 500]); hold on;
h = worldmap([g.lat(1)-2 g.lat(end)+2], [g.lon(1)-2 g.lon(end)+2]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
surfm(g.lat,g.lon,dcoast)
plotm(coastlat,coastlon,'k');
colorbar

Define season

season = nan(size(time));
season(time<datetime('1-March-2018')|time>=datetime('1-Dec-2018'))=1;
season(time>=datetime('1-March-2018')&time<datetime('1-May-2018'))=2;
season(time>=datetime('1-May-2018')&time<datetime('1-Sep-2018'))=3;
season(time>=datetime('1-Sep-2018')&time<datetime('1-Dec-2018'))=4;
Illustration
figure('position',[0 0 1000 300]); hold on;
plot(g.time,mean(gd.dens_est,'omitnan'))
plot(time,season*10)

Define night/day time

file='./data/ECMWF/2018_ssrdc.nc'; %ncdisp(file);
ssrdc = permute(flip(ncread(file,'ssrdc'),2) , [2 1 3]);
daynight = ssrdc>40000;

Construct turbine sweep area

[G,ID] = findgroups(tim.grid_id);
turbineSweptArea=nan(1,g.nlm);
turbineSweptArea(ID) = splitapply(@nansum,tim.turbine_sweptArea.*tim.Number_of_turbines,G);
turbineSweptAreaMap = nan(size(g.latlonmask));
turbineSweptAreaMap(g.latlonmask) = turbineSweptArea;
Display data
figure('position',[0 0 1400 700]);
subplot(1,2,1); hold on;
h = worldmap([g.lat(1)-2 g.lat(end)+2], [g.lon(1)-2 g.lon(end)+2]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
bordersm('countries','facecolor',[228 238 243]./255);
tmp = scatterm(tim.Latitude,tim.Longitude,tim.Number_of_turbines*2,tim.Number_of_turbines.*tim.turbine_sweptArea,'filled','MarkerFaceAlpha',.2);
c=colorbar('southoutside'); c.Label.String='Total Swept Area [m^2]';
legend(tmp,'size=Number of windturbine')
subplot(1,2,2); hold on;
h = worldmap([g.lat(1)-2 g.lat(end)+2], [g.lon(1)-2 g.lon(end)+2]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
bordersm('countries','facecolor',[228 238 243]./255);
surfm(g.lat,g.lon,turbineSweptAreaMap)
c=colorbar('southoutside'); c.Label.String='Total Swept Area [m^2]';
% colormap(brewermap([],'YlGn'));
colormap(greens); caxis([0 1000000])
Geographical characteristic
tmp = table({'France','Germany','Benelux'}',...
[sum(area(inCountry(:,:,1))) ; sum(area(inCountry(:,:,2))) ; sum(area(inCountry(:,:,3)))],...
[sum(tim.Number_of_turbines(strcmp(tim.Country,'France'))); sum(tim.Number_of_turbines(strcmp(tim.Country,'Germany'))); sum(tim.Number_of_turbines(strcmp(tim.Country,'Netherlands')))+sum(tim.Number_of_turbines(strcmp(tim.Country,'Belgium')))+sum(tim.Number_of_turbines(strcmp(tim.Country,'Luxembourg')))],...
[nansum(turbineSweptAreaMap(inCountry(:,:,1))) ; nansum(turbineSweptAreaMap(inCountry(:,:,2))); nansum(turbineSweptAreaMap(inCountry(:,:,3)))],...
[sum(tim.Total_power(strcmp(tim.Country,'France'))); sum(tim.Total_power(strcmp(tim.Country,'Germany'))); sum(tim.Total_power(strcmp(tim.Country,'Netherlands')))+sum(tim.Total_power(strcmp(tim.Country,'Belgium')))+sum(tim.Total_power(strcmp(tim.Country,'Luxembourg')))],...
'VariableNames',{'Region','Region Area','Number of windturbine','Total Swept Area','Total Power'});
tmp.("Swept Area/Region Area") = tmp.("Total Swept Area") ./ tmp.("Region Area");
tmp
tmp = 3×6 table
 RegionRegion AreaNumber of windturbineTotal Swept AreaTotal PowerSwept Area/Region Area
1'France'5.3129e+0576144.7491e+071.6036e+0789.3882
2'Germany'3.5154e+05304721.4878e+085.3427e+07423.2077
3'Benelux'6.8337e+0426401.4756e+075150260215.9290
figure('position',[0 0 1400 600]);
scatter(dcoast(:),turbineSweptAreaMap(:)./area(:),[],reshape(sum(inCountry.*reshape(1:3,1,1,[]),3),1,[]),'filled'); l=lsline; l.Color='k'; l.LineWidth=2;
ylabel('Swept Area/Region Area'); xlabel('Distance to coast [km]')
c=colorbar; c.Ticks=0:3; c.TickLabels=['Other';tmp.Region]; box on; grid on;

Bird at risk

We compute here the total number of bird at risk on the same grid (1hr, 0.25°)

Interpolate the ratio on a grid

First, we compute the ratio of bird flying at windturbine height in time and space
g_ratio = predict(compactmdl, [repmat(g.lat2D(g.latlonmask),g.nt,1) repmat(g.lon2D(g.latlonmask),g.nt,1) g.NNT(:) repelem(datenum(g.time-g.time(1)),g.nlm,1) ones(numel(g.NNT(:)),2)]);
g_ratio = reshape(g_ratio,size(g.NNT));
g_ratio(g_ratio<0)=0;
g_ratio(g_ratio>1)=1;
Same for the abscolute error
g_ratio_abs = predict(compactmdlabserror, gd.dens_est(:));

Compute birdflow

Then, we compute the bird flow [bird/m^2/hr] as
birdFlow_15min = nan(g.nlat,g.nlon,g.nt);
birdFlow_15min(repmat(g.latlonmask,1,1,g.nt)) = gd.dens_est * 1e-6 .* g_ratio .* sqrt(guv.u_est.^2 + guv.v_est.^2)*60*60;
UNIT check: bird/km^2 * km^2/m^2 * - * 1/m * m/s * s/hr -> bird/m^2 * 1/m * m/hr -> bird/m^2/hr

Downscale birdflow to 1hour

Convert time resolution of birdflow to match energy map. Average for the previous hour as wind is also the average over the previous hour.
[~,Locb]=ismember(dateshift(g.time,'end','hour'),time);
birdFlow = nan(g.nlat,g.nlon,numel(time));
for i=1:numel(time)
birdFlow(:,:,i) = nanmean(birdFlow_15min(:,:,Locb==i),3);
end
We so the same for bird flight direction
birdDir_tmp = nan(g.nlat,g.nlon,g.nt);
birdDir_tmp(repmat(g.latlonmask,1,1,g.nt)) = angle(guv.u_est + 1i*guv.v_est);
birdDir = nan(g.nlat,g.nlon,numel(time));
for i=1:numel(time)
birdDir(:,:,i) = meanangle(birdDir_tmp(:,:,Locb==i),3);
end
% histogram(abs(cos(angle(wu + 1i*wv)-birdDir)))

Birds at Risk

Compute the bird at risk [bird/hr] as bird flow [bird/m^2/hr] x total swept area [m^2].
Suming require to divide by 4 because timstep is 15min (4 in 1 hours)
birdAtRisk = birdFlow .* turbineSweptAreaMap;
disp(['Total Number of bird at risk: ' num2str(round(nansum(birdAtRisk(:)/4))/1e6) ' Millions.'])
Total Number of bird at risk: 52.1155 Millions.

Correction for wind direction map

First, the difference of wind direction and bird flight direction reduces the area with turbines base on the projection of swpt area plane perpendicular to bird direction (perp. to wind direction)
We first need to compute the windspeed u and v component seperatly
file='./data/ECMWF/2018_srf.nc'; %ncdisp(file);
w_time = datetime('2018-01-01 00:00'):1/24:datetime('2018-12-31 23:00'); % ncread(file,'time')
w_lat=flip(double(ncread(file,'latitude')));
w_lon=double(ncread(file,'longitude'));
w_u100 = permute(flip(ncread(file,'u100'),2) , [2 1 3]);
w_v100 = permute(flip(ncread(file,'v100'),2) , [2 1 3]);
w_u10 = permute(flip(ncread(file,'u10'),2) , [2 1 3]);
w_v10 = permute(flip(ncread(file,'v10'),2) , [2 1 3]);
we assume here a similar height accross space and time
height_mean = sum(tim.Hub_height.*tim.Number_of_turbines) ./ sum(tim.Number_of_turbines);
alpha = (log(w_u100)-log(w_u10)) ./ (log(100)-log(10));
wu_height = w_u100 .* (height_mean/100).^alpha;
alpha = (log(w_v100)-log(w_v10)) ./ (log(100)-log(10));
wv_height = w_v100 .* (height_mean/100).^alpha;
Fu = griddedInterpolant({w_lat,w_lon,datenum(w_time)},wu_height,'linear','nearest');
wu = Fu({g.lat,g.lon,datenum(time)});
Fv = griddedInterpolant({w_lat,w_lon,datenum(w_time)},wv_height,'linear','nearest');
wv = Fv({g.lat,g.lon,datenum(time)});
clear Fu Fv wv_height wu_height w_* alpha
Display data
figure('position',[0 0 1400 700]);
h = worldmap([g.lat(1)-2 g.lat(end)+2], [g.lon(1)-2 g.lon(end)+2]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
bordersm('countries','facecolor',[228 238 243]./255);
quiverm(g.lat2D,g.lon2D,mean(wu,3),mean(wv,3))