matlab读取GRACE数据

本科测绘专业，研究生大地测量学，但是研究生跟的老师偏地球物理，没办法只有恶补知识！！但是最近有在学习gmt，但是令人头秃。。

本来是冲着gmt强大的绘图能力的。

海洋温度
但是我在学习中发现matlab也能做到，所以我在这里记录一下。

matlab要绘制地图，首先要配置一个东西：m_map，大家可以自己百度一下配置教程。之后就可以在matlab中绘制带有经纬度信息的地理图。

m_map自带插图
现在要解决的问题是，怎样在matlab显示grcae空间分布的数据？其实我现在才发现可以看demo中的代码修改，但是我前几天看到了冯伟大佬的软件包才白嫖到一点东西。贴上代码：

A=ncread('CSR_GRACE_GRACE-FO_RL06_Mascons_all-corrections_v02.nc','lwe_thickness');%读取全球CSR重力场数据

【这个数据是在CSR官网下载的，0.25°的分辨率】

% grid_data=A(:,:,180:190);%读取2020年1月至2020年11月的全球数据

%%提取研究流域内部的重力变化数据

Volga=zeros(620,190);

YH=textread('kong.txt');【这个是研究区域的边界，需要换成自己的】

x2=1:1:720;

y2=1:1:1440;

[X2,Y2]=meshgrid(x2,y2);

a11=reshape(X2,1440*720,1);

a12=reshape(Y2,1440*720,1);

for i=1:190

M=A(:,:,i);

a13=reshape(M,1440*720,1);

G=[a12,a11,a13];

% G=trend_global1;

M4=reshape(G(:,3),1440,720);

hold on

imagesc(M4)

scatter((YH(:,2))*4+360,YH(:,1)*4)

in=inpolygon(G(:,1),G(:,2),YH(:,1)*4,(YH(:,2))*4+360);   %%0的逻辑值代表点在流域多边形的内部

EWH1=G(in,:);                                            %%找出位于流域多边形内部的点

lon=EWH1(:,1)/4;

lat=(EWH1(:,2)-360)/4;

Volga(:,i)=EWH1(:,3);                         %%将单位换成mm

disp(['已读取第',num2str(i),'层数据！']);

if(i==190)

    disp('数据读取完毕！！！！--------------------->>>>>');

end

end

lona=EWH1(:,1)/4;

lata=(EWH1(:,2)-360)/4;

grid_data=Volga(:,179:190);%%提取研究区域内部的水储量变化

[不要问为什么是这个Volga，只是一个变量]

color=max(grid_data);

title1=['1','2','3','4','5','6','7','8','9','10','11'];

%%绘制每月的变化情况

pp=textread('pp.txt');

pp1=textread('bound.txt');

for i=1:11

colorbar_value=50;

colorbar_unit='mm';

title_string=['2020年',num2str(i),'水储量'];

filename='name1';

%%全球经纬度格网

% lon = 0.125:0.25:359.875;

% lat = 89.875:-0.25:-89.875;

%%鄱阳湖其余经纬度格网化

lon = 112.5:0.25:120;

lat = 26.5:0.25:31.25;

[LON,LAT]=meshgrid(lon,lat);

% gmt_grid2map(grid_data,colorbar_value,colorbar_unit,title_string,filename)

m_proj('Equidistant Cylindrical','long',[112.5 120],'lat',[26.5 31.25]);

% m_proj('Robinson','long',[0 360],'lat',[-90 90]);

C=reshape(grid_data(:,i),31,20);

B=fliplr(C);  %%矩阵上下翻转

subplot(4,3,i)  【循环绘制子图】

m_pcolor(LON,LAT,B');

hold on

m_scatter(pp(:,1),pp(:,2),1,'b')

m_plot(pp1(:,1),pp1(:,2))

shading flat;

m_coast('linewidth',1,'color','k');

m_grid('xtick',5,'ytick',5,'tickdir','in','xlabeldir','middle',...

    'TickLength',0.008,'LineWidth',1.,'FontName', 'Helvetica','FontSize',9,'fontweight','bold');

% title([title_string],'fontsize',12,'FontName', 'Helvetica','fontweight','bold');

title(title_string)

caxis([-colorbar_value,colorbar_value]);

colormap('jet')

h=colorbar('v','FontSize',12,'fontweight','bold');

set(get(h,'title'),'string',colorbar_unit,'FontSize',12,'fontweight','bold');

end

kj=[lona,lata,grid_data];

% set(gcf, 'Color', 'white'); % white bckgr

% export_fig(gcf, ...      % figure handle

%     filename,... % name of output file without extension

%     '-painters', ...      % renderer

%     '-jpg', ...           % file format

%     '-r300' );             % resolution in dpi

【后面注释掉的代码是一个可以导出高清图片的代码】

这样得到的结果如下图：

刚刚最后一处的代码，其实是用于导出图片的。matlab系统导出的图片格式分辨率不是很清楚，而这个函数可以导出更好的图片格式。【export_fig - File Exchange - MATLAB Central (mathworks.com)】。大家可以试试这个导出的函数

 作者：我是水怪的哥 https://www.bilibili.com/read/cv11385119?spm_id_from=333.999.0.0 出处：bilibili