高精度全球地势topography的MATLAB作图（基于ETOPO1数据）

cfltwy

   
     由于该数据精度相当高，原始的Etopo数据是ETOPO1_Bed_g_gmt4.grd格式，有近400Mb,  普通电脑（8G内存以下）根本无法正常画出全球的地势contour图，所以在这里截取了包括中国在内的部分经纬度etopo1数据（北纬60-160   东经10-60）  。附件给出了截取的经纬度及etopo1数据，是MATLAB基本.mat格式。    程序中国界和省界的画法有在我第一篇帖子里讲到，可下载。

clc
clear all
close all

load etopo etopo
load lat lat
load lon lon
etopo1=(etopo(1:4:end,1:4:end))';
lat=lat(1:4:end);lon=lon(1:4:end);
% etopo1=etopo';
[xx,yy]=meshgrid(lon,lat);
  h1 = figure;
figure(h1);

load coast;
plot(long,lat,'-k','LineWidth', 1);%画世界陆地界线

  hold on

geoshow('bou1_4p.shp','facecolor',[1 1 1])    %%%%%%%%%%%%国界      boul2_4p 是省界
axis([60 160 10 60]);
% axis on
% mb=shaperead('bou2_4p.shp');% 省界
% bou2_4lx=[mb(:).X];% 提取经度
% bou2_4ly= [mb(:).Y];% 提取纬度
% m_proj('Mercator','lon',[60,150],'lat',[10,60]);% 投影
% m_plot(bou2_4lx,bou2_4ly);%绘图，查看边界是否正确，注意其实国界在省界之内，不必提取国界也可以由省界表示出国界来。
alpha(0);   %%透明度函数
% axis([60 160 10 60])
x=60.02*ones(1000,1);
y=linspace(10,60,1000);
plot(x,y,'k')
hold on
y=10.02*ones(1000,1);
x=linspace(60,160,1000);
plot(x,y,'k')
hold on
contourf(xx,yy,etopo1,'k','linestyle','none');   %%%画出等值线

axis([60 160 10 60])
%  m_proj('mercator','lon',[80 140],'lat',[10 60]);
% m_tbase('contourf');
% m_grid('linestyle','none','tickdir','out','linewidth',1.5);
%   smooth(data,10,'moving');

h=colorbar ;
xlabel(h,'m');  %定义南边彩色条水平标签
set(h,'position',[0.87  0.45  0.035  0.45]);%%%%定义彩色条位置及其长 宽 高
set (gca,'position',[0.1,0.15,0.75,0.75] );  %%%%%%%%%%固定框图长 宽 高
hold on


%  plot(122.5,52,'k*');plot(116.2,40.3,'k*');plot(109,19.4,'k*');plot(114.5,31,'k*');plot(114,29.5,'ko');plot(117,39.5,'ko');
%  imagesc(0:2.5:357.5,-87.5:2.5:90,(ncep_data(1:144,2:73,35*4))'); set(gca,'YDir','normal');
%  hold on
  plot(114.5,31,'ko','MarkerSize',6,'MarkerFaceColor','k');
  hold on
  plot(117.5,40,'ko','MarkerSize',6,'MarkerFaceColor','k');

% set(h,'ShowText','on','LevelStep',6,'TextStep',get(h,'LevelStep')*2)    %%%命令set和get是设置和获得句柄对象属性值的函数。用法是set(句柄,属性1,属性值,属性2,属性值,……)属性值=get(句柄,属性)h是句柄，这里就是指向某等高线图的句柄。你可以用get(h)看到其所有属性的属性值。showtext、textstep、levelstep是该对象的属性。ShowText为on就是显示等高线的值，off就是不显示等高线的值，默认为off。LevelStep为m就是显示值为m的倍数的等高线，显然m越大，所绘制的等高线就越稀疏。TextStep为n就是显示值为n的倍数的等高线的值
% set(h,'ShowText','on','LevelStep',6,'TextList',[300,252])
% set(h,'ShowText','on','TextList',[1,6]*1e-2)%显示1e-2和6e-2两个等高线的值colorbar


% set (gca,'position',[0.1,0.15,0.75,0.75] );  %%%%%%%%%%固定框图长 宽 高

axis([60 160 10 60])
%axis([Lon(43) Lon(57) Lat(29) Lat(15)])
hold off;

% set(gca,'XLim',[100 135],'YLim',[10 55], ...
%     'XTick',[100 105 110 115 120 125 130 135], ...
%     'Ytick',[10 15 20 25 30 35 40 45 50 55]);
xlabel('Longitude (^oE)','fontsize', 11);
ylabel('Latitude (^oN)','fontsize', 11)
           set(gcf,'color','white')    %背景底色为白色。

xyan88

thank you for your share

xyan88

警告: 未找到变量 'lat'。
> In topography at 6
警告: 未找到变量 'lon'。
> In topography at 7
未定义函数或变量 'lat'。

出错 topography (line 9)
lat=lat(1:4:end);lon=lon(1:4:end);
请问什么原因啊

cfltwy

xyan88 发表于 2016-5-14 12:26
警告: 未找到变量 'lat'。
> In topography at 6
警告: 未找到变量 'lon'。

下载的.mat文件是etopo1 lat.mat,  你自己改成对应的lat.mat文件就OK了。

cfltwy

xyan88 发表于 2016-5-14 12:18
thank you for your share

  如果有什么错误，欢迎指正！

Will

十分感谢，正需要这个地形高度文件。

Will

@cfltwy 楼主，范围貌似搞错了，应该是东经60-160，北纬10-60.

lucy1009

谢谢楼主！

lucy1009

thanks very mucj!

lucy1009

congratulate for you ! thank you very much