运用下面这期博文中的格网数据
绘制全球各大洲典型流域的时间序列图-CSDN博客https://blog.csdn.net/WH_G_Y/article/details/139972124?spm=1001.2014.3001.5501大地水准面阶方差计算公式:
计算中所用到两个转换函数:
1、“GenerateC2030C2021TransIndex”函数 功能:调换球谐系数的排列顺序(例如按阶或者按次升序排列)
function [C2021toC2030_index,C2030toC2021_index]=GenerateC2030C2021TransIndex(lmax)
% The CSR SHC storage: C20,C30,C40,...,C21,S21,C22,S22,C31,S31,C32,S32,C33,S33,...
% We define the above as "C2030" format
% Another storage(such as ITSG):C20,C21,S21,C22,S22,C30,C31,S31,C32,S32,C33,S33,....
% We define the above as "C2021" format
index=1;
lm=zeros(1,2);
for i=2:lmax
for j=0:i
if j==0
lm(index,:)=[i,j];
index=index+1;
else
lm(index:index+1,:)=[i,j;i,j];
index=index+2;
end
end
end
C2021toC2030_index=[find(lm(:,2)==0);find(lm(:,2)~=0)];
lm_=zeros(1,2);index_=1;
for i=2:lmax
for j=0:i
if j==0
lm_(index_,:)=[i,j];
index_=index_+1;
end
end
end
for i=2:lmax
for j=0:i
if j~=0
lm_(index_:index_+1,:)=[i,j;i,j];
index_=index_+2;
end
end
end
C2030toC2021_index=zeros(size(lm_,1),1);ind=1;
for ii=1:lmax
for jj=1:size(lm_,1)
if lm_(jj,1)==ii
C2030toC2021_index(ind,1)=jj;
ind=ind+1;
end
end
end
end
2、“Tran_vector ”函数 功能:将方阵存储的球谐系数展开为列向量
function [CS] = Tran_vector(cs,lmax,num,index)
CS=zeros((lmax+1)^2,num);
if ndims(cs)==2
sc=gmt_cs2sc(cs);
y_=sc2SHCvector(sc);%Do not contain C20 and degree_1 %球谐系数存放为s、c交替出现
CS(5:(lmax+1)^2,:)=y_(index,1);%将其转换由c递增转换为s递增
CS(1:2,:)=sc(1:2,lmax+1);CS(3,:)=sc(2,lmax+2);CS(4,:)=sc(2,lmax);
else
for i=1:num
sc=gmt_cs2sc(reshape(cs(i,:,:),lmax+1,lmax+1));
y_=sc2SHCvector(sc);%Do not contain C20 and degree_1 %球谐系数存放为s、c交替出现
CS(5:(lmax+1)^2,i)=y_(index,1);%将其转换由c递增转换为s递增
CS(1:2,i)=sc(1:2,lmax+1);CS(3,i)=sc(2,lmax+2);CS(4,i)=sc(2,lmax);
end
end
end
核心计算和绘图程序:
clearvars -except
addpath E:\Data\result\grid
load grid_CSR.mat;load grid300_CSR.mat;load gridP4M6_CSR.mat;load grid300P4M6_CSR.mat;
name{1,1}='_CSR';name{2,1}='300_CSR';
name{3,1}='P4M6_CSR';name{4,1}='300P4M6_CSR';
k=size(name,1); num=size(grid_CSR,3);
for jj=1:k
cs=gmt_grid2cs(eval(['grid' name{jj,1};]),60); %%转换为质量系数
for ii=1:num
temp(:,:)=cs(ii,:,:);
cs_temp(ii,:,:)=gmt_mc2gc(temp); %%质量系数转换为球谐系数
end
switch jj
case 1
cs_CSR=cs_temp;
case 2
cs300_CSR=cs_temp;
case 3
csP4M6_CSR=cs_temp;
otherwise
cs300P4M6_CSR=cs_temp;
end
end
lmax=60; [~,C2030toC2021_index]=GenerateC2030C2021TransIndex(lmax);
cs1 = Tran_vector(cs_CSR,lmax,num,C2030toC2021_index); %%将61*61方阵的球谐展开为列向量
cs2 = Tran_vector(cs300_CSR,lmax,num,C2030toC2021_index);
cs3 = Tran_vector(csP4M6_CSR,lmax,num,C2030toC2021_index);
cs4 = Tran_vector(cs300P4M6_CSR,lmax,num,C2030toC2021_index);
csm(:,1)=mean(cs1,2);csm(:,2)=mean(cs2,2); csm(:,3)=mean(cs3,2);csm(:,4)=mean(cs4,2);
for j=1:k
index=5;
for i=2:lmax
Temp(i,j)=6371e3*sqrt(mean(csm(index:index+2*i,j).^2));
index=index+2*i+1;
end
end
close all;
l=2:lmax;th=1.5;
colorS=[0 1 0;0 0 1;1 0 1;0.1 0.1 0.5];
figure(1);
for i=1:k
plot(l,abs(Temp(2:end,i)),'color',colorS(i,:),'Linewidth',th);hold on;
end
grid on;
set(gca,'yscale','log');set(gca,'Fontsize',10,'ylim',[5e-8,1e-3]);
xticks([0 10 20 30 40 50 60]);
xlabel('Degree','Fontname','Time New Roman');
ylabel('Geoid degree variance','Fontname','Time New Roman');
l1=legend('Unfiltered','G300','P4M6','G300+P4M6');
set(l1,'box','off',"FontSize",10);
结果图谱:
欢迎评论区或私信交流,多多点赞收藏,感谢支持。
参考文献:
Wahr,John.Time-variable gravity from GRACE: First results[J].Geophysical Research Letters, 2004, 31(11):293-317.DOI:10.1029/2004GL019779.