地理数据中的分辨率转换

数据分辨率问题

---

气象海洋数据在实际应用中，常常涉及到重采样，即分辨率的提高或降低等操作。本文提供了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