数据分辨率问题
气象海洋数据在实际应用中,常常涉及到重采样,即分辨率的提高或降低等操作。本文提供了matlab以及python的样例程序,以降低(网格平均)或提高(线性插值)数据的分辨率。
1. 高分辨率 ——> 低分辨率
可以使用循环逐个网格进行操作,但循环次数过多,存在效率低下的问题。
% --- 需要的分辨率0.25°, 以及经纬度网格点
deg = 0.25;
lat_era = 16:deg:47.75;
lon_era = 96:deg:127.75;
% --- 原始数据的经纬度数据
lat = ncread(full_path,'lat');
lon = ncread(full_path,'lon');
for ii=1:length(lat_era)
for jj=1:length(lon_era)
ix1 = find((lat_era(ii)-deg/2<=lat & lat<=lat_era(ii)+deg/2));
ix2 = find((lon_era(jj)-deg/2<=lon & lon<=lon_era(jj)+deg/2));
if ~isempty(ix1) & ~isempty(ix2)
t2(jj,ii) = nanmean(nanmean(T2(ix2,ix1))) - 273.15;
d2(jj,ii) = nanmean(nanmean(D2(ix2,ix1))) - 273.15;
v10(jj,ii)= nanmean(nanmean(V10(ix2,ix1)));
u10(jj,ii)= nanmean(nanmean(U10(ix2,ix1)));
end
end
end
可以基于矩阵操作, 对数据进行维度转换再平均, 由精细分辨率转换到粗分辨率.
但当原始分辨率网格不能被降低倍数整除时,需要提前进行降低分辨率之后的经纬度网格构建,对边缘数据进行填充,再进行矩阵维度的变换计算
Matlab版本
function Z=finterpHtoL(A,X)
% A -- original data
% X -- times
[nx,ny] = size(A);
X2 = X*X;
A = reshape(A,nx,X,ny/X);
A = permute(A,[2,1,3]);
A = reshape(A,X2,nx*ny/X2);
mask = ~isnan(A);
mask = sum(mask);
mask(mask<(X2/2))=nan;
A = nanmean(A);
A(isnan(mask))=nan;
Z = reshape(A,nx/X,ny/X);
return
clear;clc
A = randn(600,500);
B = finterpHtoL(A,5); % --- 分辨率降低5倍,数据维度[120,100]
figure;contourf(A');colorbar
figure;contourf(B');colorbar
Python版本
import numpy as np
def ConvertLow(data, times):
[nx, ny, nt] = data.shape
data = np.reshape(data, (int(nx/times), times, ny, nt))
data = data.transpose(0,2,1,3)
data = np.reshape(data, (int(nx/times), int(ny/times), times*times, nt))
res = np.mean(data, 2)
return res
2. 低分辨率 ——> 高分辨率
与高分数据转换为低分数据一样,我之前处理的时候都是使用matlab中的griddata函数进行双线性插值计算,效率很低.
使用python中的上采样工具,计算速度更快,可以得到同样的效果.
同样地,当获取的更精细的分辨率不是粗分辨率网格的整数倍时,python需要进行边缘网格填充处理. 使用matlab中的griddata方法不存在这一问题
Matlab版本
lat_era = 16:0.25:47.75;
lon_era = 96:0.25:127.75;
[lon_eram,lat_eram] = meshgrid(lon_era,lat_era);
lonh = 110.35:0.05:116.70;
lath = 30.8:0.05:37.15;
[long,latg] = meshgrid(lonh,lath);
for ii = 1:size(zo,3)
q(:,:,ii)= griddata(lon_eram,lat_eram,double(zo(:,:,ii)),long,latg,'linear');
end
使用interp2函数进行双线性插值比griddata函数高效得多。
%---Interp Grid
lonh = 109:0.03:117.97;
lath = 20:0.03:25.97;
[Xq, Yq] = meshgrid(lath,lonh);
%---Origin grid
lat_x = 0.0313:0.0625:64.9688;
lon_x = 60.0313:0.0625:159.9688;
[X, Y] = meshgrid(lat_x,lon_x);
%---data: the origin data
data_interp = interp2(X, Y, data, Xq, Yq, 'linear')
Python版本
from torch import nn
def ConvertHigh(data):
up = nn.Upsample(scale_factor=5, mode='bilinear', align_corners=True)
res = up(data)
return res
2023/05/06更新
"""
Funtion: This file contains the interpolate method of data. Specifically including the following content,
*** create: 2023-04-21 ***
-- convert the origin 815*430 dimension to 256*256 or 128*128,
-- convert the interpolate 256*256 or 128*128 to the origin 815*430,
Author: Muyeqingfeng
e-mail: muyeqingfeng@126.com
"""
import numpy as np
from scipy import interpolate
class Conver_Resolution(object):
def __init__(self):
super(Conver_Resolution, self).__init__()
# the origin grid
self.x_origin = np.arange(0, 815)
self.y_origin = np.arange(0, 430)
# set the interpolate grid points, convert the data dimensions to 256*256
self.x_interp = np.linspace(0, 815, 256)
self.y_interp = np.linspace(0, 430, 256)
def conver256(self, data):
"""the function to convert single 2D data, from 815*430 to 256*256"""
# define the interpolate function
interp_func = interpolate.interp2d(self.x_origin, self.y_origin, data)
# calculate the value of interpolate grid points
data_interp = interp_func(self.x_interp, self.y_interp)
return data_interp
def conver815(self, data):
"""the function to convert single 2D data, from 256*256 to 815*430"""
# define the interpolate function
interp_func = interpolate.interp2d(self.x_interp, self.y_interp, data)
# calculate the value of interpolate grid points
data_interp = interp_func(self.x_origin, self.y_origin)
return data_interp