本人:B站同步
我们知道,任何球面上的一组观测值可以球谐展开至一定的阶数,见如下的公式:
而为了将GLDAS提取的水储量与GRACE能够进行对比,需要进行同样的截断滤波处理。本专栏以GLDAS数据的土壤水为例,分析如何进行球谐展开并滤波。主要代码如下:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GLDAS expand
% Chistrong Wen
% 2022-12-27
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GLDAS INFO
file = 'GLDAS_VIC10_M.A201710.021.nc4';
ncdisp(file)
%% READ SOIL MOISTURE
lon = ncread(file,'lon');
lat = ncread(file,'lat');
[lon,lat] = meshgrid(lon,lat);
%%
lon = -179.5:1:179.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);
%% SoilMoi0_30cm_inst SoilMoi_depth2_inst SoilMoi_depth3_inst
soil1 = ncread(file,'SoilMoi0_30cm_inst');
soil2 = ncread(file,'SoilMoi_depth2_inst');
soil3 = ncread(file,'SoilMoi_depth3_inst');
soil_sum = soil1+soil2+soil3;
out_soil.lon = lon;
out_soil.lat = lat;
fill = nan(30,360);
out_soil_temp = [fill;soil_sum'];
out_soil.rg = out_soil_temp;
figure
wzq_plot(out_soil)
title('raw data')
grid = out_soil.rg;
degree_max = 60;
[cs] = gmt_grid2cs(grid,degree_max);
destrip_method = 'SWENSON';
radius_filter = 300;
type = 1;
[grid_filter]=gmt_cs2grid(cs,radius_filter,type,destrip_method)
lon = -179.5:1:179.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);
figure
out_soil_filter.lon = lon;
out_soil_filter.lat = lat;
out_soil_filter.rg = grid_filter;
wzq_plot(out_soil_filter)
title('degree 60')
%% truncted to 96
degree_max = 96;
[cs] = gmt_grid2cs(grid,degree_max);
destrip_method = 'SWENSON';
radius_filter = 300;
type = 1;
[grid_filter]=gmt_cs2grid(cs,radius_filter,type,destrip_method)
lon = -179.5:1:179.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);
figure
out_soil_filter.lon = lon;
out_soil_filter.lat = lat;
out_soil_filter.rg = grid_filter;
wzq_plot(out_soil_filter)
title('degree 96')
%% smooth at 500 km
degree_max = 60;
[cs] = gmt_grid2cs(grid,degree_max);
destrip_method = 'SWENSON';
radius_filter = 500;
type = 1;
[grid_filter]=gmt_cs2grid(cs,radius_filter,type,destrip_method)
lon = -179.5:1:179.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);
figure
out_soil_filter.lon = lon;
out_soil_filter.lat = lat;
out_soil_filter.rg = grid_filter;
wzq_plot(out_soil_filter)
title('Gaussian 500 km')
下面是截断滤波后的结果:我们看到截断滤波后,信号强度有所衰减,而随着截断阶数越大,信号越接近真实信号强度。而高斯平滑的半径越大,信号也衰减。
欢迎大家收藏点赞❤❤❤