cubic interp2 on latitude & longitude data. wind comparison

4 次查看(过去 30 天)
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')

采纳的回答

Mathieu NOE
Mathieu NOE 2023-7-17
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--');
  2 个评论
Tiziano Bagnasco
Tiziano Bagnasco 2023-7-21
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 :)

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Surface and Mesh Plots 的更多信息

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by