请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 27509|回复: 13

[程序设计] 用ETOPO2数据,通过m_map画水深分布图

[复制链接]

新浪微博达人勋

发表于 2017-5-5 23:54:04 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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);%将坐标从-18090,变为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世界地图底图

m_map世界地图底图
                 图1   画出worldmap



我使用的数据:ETOPO2v2g_f4.nc,下载地址为
这个错误是因为,没有sort将经度按从小到大的顺序排列的结果(图2)。
      

画错的图

画错的图
              图2   经度未排序时画出的worldmap

下面是最终画出的水深分布图(图3),和某一区域的水深图(图3)。
在程序中可以global_water_depth1.m 中,可以选择区域,比如(图4)。
      

水深分布图

水深分布图
                    图3  最终画出的水深分布图
      

某一区域示例

某一区域示例
                        图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')

global_water_depth1.m

2.36 KB, 下载次数: 54, 下载积分: 金钱 -5

评分

参与人数 1金钱 +10 贡献 +5 收起 理由
二爷名声在外 + 10 + 5 很给力!

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-5-6 09:35:04 | 显示全部楼层
colorbar不错~~
不过180E那似乎还是有条线
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-5-7 11:17:30 | 显示全部楼层
ljchen1989 发表于 2017-5-6 09:35
colorbar不错~~
不过180E那似乎还是有条线

好像俄罗斯的楚科奇半岛那是有条线;colobar是GEBCO地形数据发布的图上取出来的颜色。http://www.gebco.net/data_and_pr ... ps/gebco_world_map/
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-12-24 22:05:39 | 显示全部楼层
要研究下这个
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-3-9 19:20:49 | 显示全部楼层
233333666666666
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-6-27 20:45:18 来自手机 | 显示全部楼层
不错,研究研究
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-3-11 11:02:05 | 显示全部楼层
不错不错,m_map强大
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-8-11 11:28:12 | 显示全部楼层
我想问下怎么我现在想下载那个数据 无法下载(自动回复:请不要使用迅雷等下载工具,点我查看下载帮助
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-4-9 16:56:07 | 显示全部楼层
本帖最后由 boyfriendxin 于 2020-4-9 16:59 编辑

请问您的colormap从那个网站怎么取的呢,怎么就知道具体的数值的呢?
%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;
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-4-9 16:59:23 | 显示全部楼层
请问怎么从那个网站获取colorbar的rgb矩阵的呢???
%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;
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表