本科测绘专业,研究生大地测量学,但是研究生跟的老师偏地球物理,没办法只有恶补知识!!但是最近有在学习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