登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 葫芦爷儿 于 2017-5-6 00:06 编辑
M_map画全球地图
遇到的问题:为了突出西太平洋的海沟分布位置,我想作出太平洋在地图的中心的水深分布图。但在m_map 画全球图时,是以0˚为中心(图1)。 通过尝试,我发现m_map的long可以写成0-360的。 但是需要注意的是,要把水深数据对应的经纬度改成0-360。 m_proj('Equidistant Cylindrical','long',[90 450],'lat',[-90 90]) %90˚E为地图的左边界 x_sss=x_ss; px=find(x_ss<90); x_sss(px)=360+x_ss(px);%将坐标从-180:90,变为180:450 x_sss=sort(x_sss); [lon,lat]=meshgrid(x_sss,y_ss); t=dep(:,1:px(end)); tt=dep(:,px(end)+1:end); dep1=[tt,t];%重新排列水深矩阵
m_map世界地图底图
我使用的数据:ETOPO2v2g_f4.nc,下载地址为 这个错误是因为,没有sort将经度按从小到大的顺序排列的结果(图2)。
画错的图
下面是最终画出的水深分布图(图3),和某一区域的水深图(图3)。 在程序中可以global_water_depth1.m 中,可以选择区域,比如(图4)。
水深分布图
某一区域示例
% 这个程序是用来画某一个研究区域水深地形分布图。 % 由于,m_map默认的是欧洲地图中心,我需要画出太平洋在地图的中央的水深图。 % 通过在气象家园上的查找资料,完成修改。 % 感谢气象家园,以及二爷名声在外 翻译的m_map使用手册。 clear; clc; clf; %读取数据 cd'/Users/cyao/Desktop/etopo2/' x=ncread('/Users/cyao/Desktop/etopo2/etopo2/ETOPO2v2g_f4.nc','x'); y=ncread('/Users/cyao/Desktop/etopo2/ETOPO2v2g_f4.nc','y'); z=ncread('/Users/cyao/Desktop/etopo2/ETOPO2v2g_f4.nc','z'); %确定画图区域 =find(y==-90); [i1]=find(y==90); [j]=find(x==-180); [j1]=find(x==180); z1=z(j:j1,i:i1); x_ss=x(j:j1)'; y1=y(i:i1); y_ss=rot90(rot90(y1)); z1=rot90(z1); dep=z1; dep(find(z1>0))=0; %由于数据太多,电脑跑不动,在降低分辨率之后,进行画图。 dep=dep(1:10:end,1:10:end); x_ss=x_ss(1:10:end); y_ss=y_ss(1:10:end); %从90?E切开 x_sss=x_ss; px=find(x_ss<90); x_sss(px)=360+x_ss(px); x_sss=sort(x_sss); [lat,lon]=meshgrid(x_sss,y_ss); t=dep(:,1:px(end));%前半段取出来 tt=dep(:,px(end)+1:end);%后半段取出来 dep1=[tt,t];%合起来 dep2=[tt,t]; %放大colorbar p_200=find(dep1>=-200);%确定200米以浅的数据位置 p_500=find(-200>dep1&dep1>=-500);%确定200米-500米的数据位置 p_1000=find(-500>dep1&dep1>=-1000);%确定500米-1000米的数据位置 p_other=find(dep1<-1000);%确定1000米后的数据位置 %放大 dep1(p_200)=dep1(p_200)*5; dep1(p_500)=(dep1(p_500)*5+1000)*10/15-1000;%通过上次放大后,在减去上一次放大的极限值1000, %根据现在长度进行局部放大,之后在加上上一段水深的 dep1(p_1000)=(((dep1(p_1000)*5+1000)*10/15-1000)+2000)*1000*3/5000-2000; dep1(p_other)=dep1(p_other)-2000; %mycoloabr my=[196 226214;159 211 206;122 199 213;83 189 209; 60 161 192;50 123 171;43 87 145;36 67111;31 43 73;14 19 36]/255; my=flipud(my); %画图 m_proj('EquidistantCylindrical','long',[90 450],'lat',[-90 90]); h=m_pcolor(lat,lon,dep1); m_gshhs_h('color',[.7.7 .7]); m_grid('linewi',2,'tickdir','on','fontsize',15); m_coast('patch',[.7.7 .7],'edgecolor','k'); hold on; %微调图像 set(h,'Linestyle','none'); colormap(my) h1=colorbar('southoutside','Direction','reverse','Ticks',[-10000,-9000,-8000,-7000,-6000,-5000,-4000,-3000,-2000,-1000],... 'TickLabels',{'8000','7000','6000','5000','4000','3000','2000','1000','500','200',},'FontSize',13) caxis([-10000,0]) h1.Label.String ={ '水深/m ','Bathymetry'}; h1.Label.FontSize=18 %保存 print('-dpng','-r600','colorbar') |