cubic interp2 on latitude & longitude data. wind comparison

Hello everyone!

I am carrying out a quantile-quantile analysis and bias correction between two horizontal wind velocity fields for the whole year 2018: ERA5 wind data (81x109x8760,hourly data) and IFREMER wind data (82x110x1460,every 6 hours) . My domain covers the majority of the South China Sea and the two datasets are given as 3D matrices that have same resolutions ( 0.25 degrees x 0.25 degrees) but different coordinates. After downsampling the ERA5 field, in order to have data every 6 hours, I want to interpolate using interp2 and therefore find the IFREMER horizontal wind speed values at the ERA5 locations. The problem is that, after downlaoding the data from the web, the two fields are not defined in the same way: i had to flip the second dimension of the ERA5 wind field in order to obtain a comparable field ( see below) . After doing so, and after interpolating, I transform both 3D fields in column vectors so that they can be represented in a scatter plot and qq-plot.

However, what i can see from the plot is that ERA5 winds ( reanalysis winds) tend to be higher than IFREMER winds ( more accurate wind coming from scatterometers). I was expecting the opposite behaviour, as demonstrated by other resarchers.

Can you spot any mistakes or calculations that were not computed properly?

Below you can see the steps that I have taken.

Cheers and thank you a lot in advance.

% Attached you can find the ERA5 wind field --> u10_ERA5_6hours (81x109x48)
%                           IFREMER wind field ---> u10_ifremer (82x110x48)
% the entire fields were too heavy to be uploaded. These are related to 48
% hours but they should be enough to understand how the code works.

% Coordinates
lon_ERA5 = [104:0.25:124]';
lat_ERA5 = [29:0.25:2]'; % descending order
lon_IFREMER = [103.8125:0.25:124.0625]';
lat_IFREMER = [1.8125:0.25:29.0625]'; %ascending order


% assign NaN to land points in IFREMER data
mask = ~isfinite(u10_ifremer) | abs(u10_ifremer) > 1e30*eps;
u10_ifremer(mask)=NaN;

% plot the wind fields for time 1 for both IFREMER and ERA5 to understand their orientations
figure;
contourf(u10_ifremer(:,:,1))
figure;
contourf(u10_ERA5_6hours(:,:,1))

% let's flip the ERA5 matrix in order to have the same orientation
u10_ERA5_6hours = flip(u10_ERA5_6hours,2); % flip the latitude
% plot for checking it
figure;
contourf(u10_ERA5_6hours(:,:,1)) % OK !!!
% let's flip ERA5 latitude
lat_ERA5_flip = flip(lat_ERA5);

% MESHGRIDS GENERATION AND INTERP2
[X,Y] = meshgrid(lat_IFREMER, lon_IFREMER);
[Xp,Yp] = meshgrid(lat_ERA5_flip,lon_ERA5);
for t=1:1460
u10_ifremer_int(:,:,t)= interp2(X,Y,u10_ifremer(:,:,t),Xp,Yp,'cubic');
end
% Reshape command
windField1 = reshape(u10_ifremer_int, [],1);
windField2 = reshape(u10_ERA5_6hours, [],1);
% Scatter plot and Q-Q plot
figure;
scatter(windField1,windField2,".") 
axis equal
pbaspect([ 1 1 1])
hold on
qqplot(windField1,windField2)
xlabel('IFREMER u10 wind component [m/s]')
ylabel('ERA5 u10 wind component [m/s]')
title('scatter plot and Q-Q plot year 2018')

Answer 1

[Mathieu NOE]

hello

I am not an expert for atmospheric data post processing, buth here are my findings

1 / seems to me flipping the data is not appropriate

% let's flip the ERA5 matrix in order to have the same orientation

% u10_ERA5_6hours = flip(u10_ERA5_6hours,2); % flip the latitude

% plot for checking it

the plot below is done without the flip operatin, otherwise it seems to me it's not oriented the same way

2/ then the scatter plot looks better

3/ I added a hist3 like (density plot) as a further refinement of the scatter plot. It tells you how dense are the data organized along the diagonal (in red). As we can see the two data sets are quite in good match

code a bit modified :

load('ifremer.mat')
load('era5_6hours_reduced.mat')
% Attached you can find the ERA5 wind field --> u10_ERA5_6hours (81x109x48)
%                           IFREMER wind field ---> u10_ifremer (82x110x48)
% the entire fields were too heavy to be uploaded. These are related to 48
% hours but they should be enough to understand how the code works.
% Coordinates
lon_ERA5 = [104:0.25:124]';
lat_ERA5 = [29:-0.25:2]'; % descending order
lon_IFREMER = [103.8125:0.25:124.0625]';
lat_IFREMER = [1.8125:0.25:29.0625]'; %ascending order
% assign NaN to land points in IFREMER data
mask = ~isfinite(u10_ifremer) | abs(u10_ifremer) > 1e30*eps;
u10_ifremer(mask)=NaN;
% plot the wind fields for time 1 for both IFREMER and ERA5 to understand their orientations
figure;
subplot(2,1,1),contourf(u10_ifremer(:,:,1))
title('IFREMER data');
caxis([-14 4])
colorbar('vert');
% let's flip the ERA5 matrix in order to have the same orientation
% u10_ERA5_6hours = flip(u10_ERA5_6hours,2); % flip the latitude
% plot for checking it
subplot(2,1,2),contourf(u10_ERA5_6hours(:,:,1)) % OK !!!
title('ERA5 data');
caxis([-14 4])
colorbar('vert');

% let's flip ERA5 latitude
lat_ERA5_flip = flip(lat_ERA5);

% MESHGRIDS GENERATION AND INTERP2
[X,Y] = meshgrid(lat_IFREMER, lon_IFREMER);
[Xp,Yp] = meshgrid(lat_ERA5_flip,lon_ERA5);

for t=1:48%1460
    u10_ifremer_int(:,:,t)= interp2(X,Y,u10_ifremer(:,:,t),Xp,Yp,'cubic');
end
% Reshape command
windField1 = reshape(u10_ifremer_int, [],1);
windField2 = reshape(u10_ERA5_6hours, [],1);
% Scatter plot and Q-Q plot
figure;
scatter(windField1,windField2,".") 
axis equal
pbaspect([ 1 1 1])


% hold on
% qqplot(windField1,windField2)
% xlabel('IFREMER u10 wind component [m/s]')
% ylabel('ERA5 u10 wind component [m/s]')
% title('scatter plot and Q-Q plot year 2018')


%% bin the data (Hist3 clone)
x = windField1;
y = windField2;
nBins = 25; %number of bins (there will be nBins + 1 edges)
XEdge = linspace(min(x),max(x),nBins+1); 
YEdge = linspace(min(y),max(y),nBins+1); 
dx = mean(diff(XEdge));
dy = mean(diff(YEdge));
Xcenter = XEdge(1:end-1)+dx/2;
Ycenter = YEdge(1:end-1)+dy/2;

[~, xBin] = histc(x, XEdge);
[~, yBin] = histc(y, YEdge);
% count number of elements per (x,y) pair
[xIdx, yIdx] = meshgrid(1:nBins, 1:nBins); 
xyPairs = [xIdx(:), yIdx(:)]; 
Z = zeros(size(xyPairs,1),1); 
for i = 1:size(xyPairs,1)
    Z(i) = sum(ismember([xBin, yBin], xyPairs(i,:), 'rows')); 
end
% Reshape nPerBin to grid size
Z = reshape(Z, [nBins, nBins]);

% plot density map 
figure
imagesc(Xcenter,Ycenter,Z);
set(gca,'YDir','normal');
xlabel('windField1');
ylabel('windField2');
axis square
xlim([-18 10]);
ylim([-18 10]);
hold on
plot((-18:10),(-18:10),'r--');

[Tiziano Bagnasco]

Thank you a lot for your answer!

However, I just double checked and the 3D era5 field attached here was flipped already! My bad!

Starting from the raw data i obtain this plot here:

Anyway, thank you for the suggestions in your answer! They helped me into refining my code. Cheers :)

[Mathieu NOE]

ok

again, my pleasure and good luck for the future