matlab论文图一的地形区域图的球形展示Version_1

matlab论文图一的地形区域图的球形展示Version_1
图片

此图来源于：

![Jieqiong Zhou, Ziyin Wu, Dineng Zhao, Weibing Guan, Chao Zhu, Burg Flemming,
Giant sand waves on the Taiwan Banks, southern Taiwan Strait: Distribution, morphometric relationships, and hydrologic influence factors in a tide-dominated environment,
Marine Geology,
Volume 427,
2020,
106238,
ISSN 0025-3227,
https://doi.org/10.1016/j.margeo.2020.106238.](https://i-blog.csdnimg.cn/direct/1af67a3eb5264436847e4fc94a92ffc0.png)

这个图的地形数据很精细，因为我画的图没展示这么精细。
底图海图可以画，然后左上角放个球形地图：
写成函数可以进行调用：

add_sphere1(cmap)

cmap指的是colormap；

调用格式： 先设置位置，在添加colormap，在调整视角即可。

%% %% add sphere  set position
axes('position',[0.02 0.49 0.42 0.42]) 
% add colormap
cmap=load('MPL_terrain.txt');
add_sphere1(cmap)
view(40,25);% view angles
% view(x,y);% x 控制左右旋转，y控制上下旋转。

结果展示：

图片

图片

图片

图片

代码：

底图海图

.rtcContent { padding: 30px; } .lineNode {font-size: 12pt; font-family: "Times New Roman", Menlo, Monaco, Consolas, "Courier New", monospace; font-style: normal; font-weight: normal; }
clc;clear;close all
% read data
file='F:\data\etopo\etopo1.nc';
lon=double(ncread(file,'x'));
lat=double(ncread( file,'y'));
h=double(ncread(file,'z'));
% 选定区域南海%
area  =[116 121 22 25];
ln =find(lon>=area(1)&lon<=area(2));
la=find(lat>=area(3)&lat<=area(4));
lon = lon(ln);
lat  = lat(la);
H = h(ln,la);
[x,y]=meshgrid(lon,lat);
x=x'; y=y';
%% %%m_pcolor画出了区域的等深线图
close all
figure;
set(0,'defaultfigurecolor','w')
set(gcf,'position',[50 50 1200 900])
m_proj('Miller','lon',[area(1) area(2)],'lat',[area(3) area(4)]);
caxis([-2500 0])% 这些必须放在m_shaderelief的前面，不然不可用
colorbar
cmap=load('MPL_terrain.txt');% add colormap
m_shadedrelief(lon,lat,H',cmap);
m_gshhs_c('patch',[0.8 0.8 0.8]);
m_grid('linest','none','xtick',[116:1:121],'ytick',[22:25],'tickdir','in',...'FontName','Times new roman','FontSize',12,'color','k','box','on');%box on off and fancy;
m_shadedrelief.m.rtcContent { padding: 30px; } .lineNode {font-size: 12pt; font-family: "Times New Roman", Menlo, Monaco, Consolas, "Courier New", monospace; font-style: normal; font-weight: normal; }
function  [Truecol,x,y]=m_shadedrelief(x,y,Z,cmap,varargin)
% M_SHADEDRELIEF Shaded relief topography in an image
%  M_SHADEDRELIEF(X,Y,Z) presents a shaded relief topography, as would
%  be seen if a 3D model was artifically lit from the side. Slopes
%  facing the light are lightened, and slopes facing away from the
%  light are darkened. X and Y are horizontal and vertical coordinate
%  VECTORS, and these should be in the same units as the height Z
%  (e.g., all in meters), otherwise the slope angle calculations will be
%  in error.
%
%  Usage notes:
%
%  (1) M_SHADEDRELIEF is a replacement for a low-level call to IMAGE
%  displaying a true-colour image so it MUST be preceded by COLORMAP and
%  CAXIS calls.
%
%  (2) Gradients have to be calculated, so Z should be a data-type accepted
%   by Matlab's GRADIENT function (currently double and single).
%   If Z is (e.g.) 'int8', or some other data-type you must use
%   M_SHADED_RELIEF(X,Y,double(Z))
%
%  (3) M_SHADEDRELIEF probably is most useful as a backdrop to maps with
%  a rectangular outline box - either a cylindrical projection, or some
%  other projection with M_PROJ(...'rectbox','on').
%
%  (4) Finally, the simplest way of not running into problems:
%     - if your elevation data is in LAT/LON coords (i.e. in a matrix where
%       each row has points with the same latitude, and each column has points
%       with the same longitude), use
%                  M_PROJ('equidistant cylindrical',...)
%     - if your elevation data is in UTM coords (meters E/N), i.e. in a matrix
%       where each row has the same UTM northing and the each column has the
%       same UTM easting, use
%                  M_PROJ('utm',....)
%
%  M_SHADEDRELIEF(...,'parameter',value) lets you set various properties.
%  These are:
%       'coords' : Coordinates of X/Y/Z:
%                     'geog' for lat/lon, Z meters,  (default)
%                     'map'  for X/Y map coordinates, Z meters
%                     'Z' if X/Y/Z are all in same units (e.g., meters)
%       'lightangle' : true direction (degrees) of light source (default
%                      -45, i.e. from the north-west)
%       'gradient': Shading effects increase with slope angle
%                   until slopes reach this value (in degrees), and are
%                   held constant for higher slopes (default 10). Reduce
%                   for smoother surfaces.
%       'clipval' : Fractional change in shading for slopes>='gradient'.
%                  0 means no change, 1 means saturation to white or black
%                  if slope is facing directly towards or away from light
%                  source, (default 0.9).
%       'nancol'  : RGB colour of NaN values (default [1 1 1]);
%       'lakecol' : RGB colour of lakes (flat sections) (default NaN)
%                   If set to NaN lakes are ignored.
%
%   IM=M_SHADEDRELIEF(...) returns a handle to the image.
%
%   [SR,X,Y]=m_SHADEDRELIEF(...) does not create an image but only returns
%   the  true-color matrix SR of the shaded relief, as well as the X/Y
%   vectors needed to display it.
%
%   Example:
%           load topo
%           subplot(2,1,1);  % Example without it
%           imagesc(topolonlim,topolatlim,topo);
%           caxis([-5000 5000]);
%           colormap([m_colmap('water',64);m_colmap('gland',64)]);
%           set(gca,'ydir','normal');
%
%           subplot(2,1,2);  % Example with it
%           caxis([-5000 5000]);
%           colormap([m_colmap('water',64);m_colmap('gland',64)]);
%           m_shadedrelief(topolonlim,topolatlim,topo,'gradient',5e2,'coord','Z');
%           axis tight
%
% Rich Pawlowicz (rich@eoas.ubc.ca) Dec/2017
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%
% Changes:
% Jan/2018 - changed outputs for flexibility
%            and added 'map' coordinate handling
% Mar/2019 - added alphamapping for out-of-map areas, started using colormap
%            local to AXES not to FIGURE.
% Apr/2019 - some parts relied on the ones-expansion; went back to meshgrid
%            for compatibility with older matlab versions (thanks P. Grahn)
% Oct/2020 - added info about how input Z must be a double
global MAP_PROJECTION MAP_VAR_LIST
lighthead=-45;
gradfac=10;
clipval=.9;
nancol=[1 1 1];
lakecol=NaN; %[.7 .9 1];
geocoords='geog';
scfac=6400000;  % Used for geo coordinates if on sphere radius 1
while ~isempty(varargin)switch lower(varargin{1}(1:3))case 'coo'switch lower(varargin{2}(1))case 'g'geocoords='geog';case 'm'geocoords='map';case {'z','u'}geocoords='z';otherwiseerror('Unknown coordinate specification');endcase 'lig'lighthead=varargin{2};case 'gra'gradfac=varargin{2};case 'cli'clipval=varargin{2};case 'nan'nancol=varargin{2};case 'lak'lakecol=varargin{2};otherwiseerror(['m_shadedrelief: Unknown option: ' varargin{1}]);endvarargin(1:2)=[];
end
% All kinds of issues dealing with coords:
% First, we need VECTOR x/y as an input to gradient function.
if isvector(x) && size(x,2)==1x=x';
elseif ~isvector(x)       % Can't be a matrixerror('Input X must be a VECTOR');
end
if isvector(y) && size(y,1)==1y=y';
elseif ~isvector(y)error('Input Y must be a VECTOR');
end
%  Now handle the case if we just give a starting and an ending point
if length(x)==2x=linspace(0,1,size(Z,2))*diff(x)+x(1);
end
if length(y)==2y=linspace(0,1,size(Z,1))'*diff(y)+y(1);
end
if strcmp(geocoords,'geog')                           %  If its Lat/Long points% Have to have initialized a map firstif isempty(MAP_PROJECTION)disp('No Map Projection initialized - call M_PROJ first!');return;end% Convert to X/Y[X,Y]=m_ll2xy(x(1,:),repmat(mean(y(:,1)),1,size(x,2)),'clip','off');[X2,Y2]=m_ll2xy(repmat(mean(x(1,:)),size(y,1),1),y(:,1),'clip','off');x=X;y=Y2;
end
% Note - 'image' spaces points evenly, so we should just check that they
% are even otherwise the image won't line up with coastlines...
if max(abs( x - linspace(x(1),x(end),length(x)) ) )/abs(x(end)-x(1)) >.005warning(['********** Image will be distorted in X direction!! use M_IMAGE to re-map? *************']);
end
if max(abs( y - linspace(y(1),y(end),length(y))' ) )/abs(y(end)-y(1)) >.005warning(['********** Image will be distorted in Y direction!! use M_IMAGE to re-map? *************']);
end
% Convert colours to uint8s
if all(nancol<=1)nancol=uint8(nancol*255);
end
% Convert colours to uint8s
if all(lakecol<=1)lakecol=uint8(lakecol*255);
end
% Get caxis
clims=caxis;
if all(clims==[0 1])   % Not setclims=[min(Z(:)) max(Z(:))];
end
% Get colormap for the current axes
% cmap1=load('GMT_drywet1.mat');
% cmap1 =(cmap1.raw_new);
% cc=colormap((cmap1));
if isempty(cmap)cmap=load('MPL_terrain.txt');
elsecmap = cmap;
end
cc=colormap((cmap));
% cc=colormap(gca);
cc2=round(cc*255);  % we need these in 0-255 range to get Truecolor
lcc=size(cc,1);
%inan=isnan(Z);
% Get slopes
% If we are using a normal ellipsoid we need to rescale
% x/y to get true slope angles
if ~isfloat(Z)warning('Your input Z matrix must be floating point');
end
if  (strcmp(geocoords,'map') || strcmp(geocoords,'geog')) && strcmp(MAP_VAR_LIST.ellipsoid,'normal')scfac=6370997;[Fx,Fy]=gradient( Z, x*scfac, y*scfac);
else[Fx,Fy]=gradient( Z, x, y);
end
% Find NaN
[inan,jnan]=find(isnan(Z) | isnan(Fx) | isnan(Fy) );
% Probable lakes
[islake,jlake]=find(Fx==0 & Fy==0);
% Convert z levels into a colormap index.
% Some iteration to discover the exact formula that matlab uses for mapping to
% color indices (from 1 to lcc)
%idx=min( floor( min(max( (Z-clims(1))/(clims(2)-clims(1)),0) ,1 )*lcc )+1,lcc);
idx=max(min( floor(   (Z-clims(1))/(clims(2)-clims(1))*lcc  )+1  ,lcc),1);
% The slope angle relative to the light direction in degrees.
Fnw=atand(imag(-(Fx+i*Fy)*exp(i*lighthead*pi/180)));
%Put an upper and lower limit on the angles
%%Fnw=min(clipval,max(-clipval,Fnw/gradfac));
Fnw=clipval*tanh(Fnw/gradfac);
% Now get the colormap for each pixel and scale the RGB value 'c'.
% If the correction is -0.1 then scale  c*(1 - |-0.1|)
% If the correction is  0.1 then scale  c*(1 - |+0.1|) + 0.1*255
%depending on the slope.
%Truecol=uint8(max(0,min(255,   reshape([cc2(idx,:)],[size(idx) 3]).*repmat(1-abs(Fnw),1,1,3)+repmat(255*Fnw.*(Fnw>0),1,1,3) ) ));
Truecol=uint8( reshape([cc2(idx,:)],[size(idx) 3]).*repmat(1-abs(Fnw),1,1,3)+repmat(255*Fnw.*(Fnw>0),1,1,3) ) ;
%Truecol=uint8(max(0,min(255,   reshape([cc2(idx,:)],[size(idx) 3]).*repmat(1+Fnw/gradfac,1,1,3) ) ));
% Colour Lakes
if any(islake) && isfinite(lakecol(1))Truecol(sub2ind(size(Truecol),islake,jlake,  ones(size(islake))))=lakecol(1);Truecol(sub2ind(size(Truecol),islake,jlake,1+ones(size(islake))))=lakecol(2);Truecol(sub2ind(size(Truecol),islake,jlake,2+ones(size(islake))))=lakecol(3);
end
% Colour the NaNs
if any(inan)Truecol(sub2ind(size(Truecol),inan,jnan,  ones(size(inan))))=nancol(1);Truecol(sub2ind(size(Truecol),inan,jnan,1+ones(size(inan))))=nancol(2);Truecol(sub2ind(size(Truecol),inan,jnan,2+ones(size(inan))))=nancol(3);
end
if strcmp(geocoords,'map')   % Have to "make invisible" the points outside the map limits.[xm,ym]=meshgrid(x,y);[HLG,HLT]=m_xy2ll(xm,ym);% Find pixels outside the limits of the actual map (if the boundary% isn't a rectangle)if strcmp(MAP_VAR_LIST.rectbox,'off')[I,J]=find(HLT<MAP_VAR_LIST.lats(1) | HLT>MAP_VAR_LIST.lats(2) | HLG<MAP_VAR_LIST.longs(1) | HLG>MAP_VAR_LIST.longs(2));elseif strcmp(MAP_VAR_LIST.rectbox,'circle')R=(xm.^2 +ym.^2);[I,J]=find(R>MAP_VAR_LIST.rhomax.^2);elseI=[];J=[];endbackcolor=uint8(get(gcf,'color')*255);if any(I)                          % if some pixels are outside the map areafor k=1:3IJ=sub2ind(size(Truecol),I,J,repmat(k,size(I)));Truecol(IJ)=backcolor(k);           % Set them to the background colourendend
elseI=[];J=[];
end
if nargout<=1if any(I)  % make pixels outside the map area transparent, if needed.alphadata=  ones(size(Truecol,1),size(Truecol,2),'logical');IJ=sub2ind(size(alphadata),I,J);alphadata(IJ)=0;Truecol=image('xdata',x,'ydata',y,'cdata',Truecol,'alphadata',alphadata,'tag','m_shadedrelief');elseTruecol=image('xdata',x,'ydata',y,'cdata',Truecol,'tag','m_shadedrelief');end
end
add_sphere1
.rtcContent { padding: 30px; } .lineNode {font-size: 12pt; font-family: "Times New Roman", Menlo, Monaco, Consolas, "Courier New", monospace; font-style: normal; font-weight: normal; }
function add_sphere1(cmap)
load topo topo topomap1    % load data
x = 0:359;                                % longitude
y = -89:90;                               % latitude
[X1,Y1]=meshgrid(x,y);
x1 = 0:0.1:360;
y1=-89:0.1:90;
[X,Y]=meshgrid(x1,y1);
topo_new = griddata(X1,Y1,topo,X,Y);
% figure;
% set(0,'defaultfigurecolor','w')
% set(gcf,'position',[50 50 1200 900])
[x,y,z] = sphere(100);          % create a sphere
s = surface(x,y,z);            % plot spherical surface
s.FaceColor = 'texturemap';    % use texture mapping
s.CData = topo_new;                % set color data to topographic data
s.EdgeColor = 'none';          % remove edges
s.FaceLighting = 'gouraud';    % preferred lighting for curved surfaces
s.SpecularStrength = 0.4;      % change the strength of the reflected light
if isempty(cmap)cmap=load('MPL_terrain.txt');
elsecmap = cmap;
end
colormap(cmap)
% light('Position',[1 0 1])     % add a light
axis square off                % set axis to square and remove axis
view([30,10])                 % set the viewing angle