该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

像这样的 ，但我想换个区域，不知道要怎么改数据

clear

clc

load('G:\weixingshuju\topo提取\kkk.mat'); %文件位置

a1 = reshape(a(:,1),4321,2161);

a2 = reshape(a(:,2),4321,2161);

a3 = reshape(a(:,3),4321,2161);

A = a1(1272:1464,805:1033); %(经度范围*12+1，纬度范围*12+1)数据像元为12*12

B = a2(1272:1464,805:1033);

C = a3(1272:1464,805:1033);

S = C';

contour(A',B',-S,[0,0],'k');

hold on

set(gca,'XTick',100:5:130);

set(gca,'Xticklabel',{'100E','105E','110E','115E','120E','125E','130E'});

set(gca,'yTick',0:5:25);

set(gca,'yticklabel',{'0N','5N','10N','15N','20N','25N'});

axis([106 122 4 23]);

xlabel(['经度/^oE']); ylabel(['纬度/^oN']); %画出地形图

load('G:\weixingshuju\weixing\Globe_SSHA_dt_ref_global_merged_msla_h_qd\data\001 (4).mat') %海面高度数据

high=ssha(377:453,425:490); %(纬度范围*4+1，经度范围*4+1)数据像元为4*4

land=find(high>1000);

high(land)=NaN;

lon1=(377:453);

lat1=lat(425:490);

[C,h]=contourf(lat1,lon1,high,'linestyle','none');

a2=load('G:\weixingshuju\topo提取\kkk.mat');

b2=reshape(a(:,3),4321,2161);

c2=b2';

d2=c2(805:1033,1272:1464);

for p=1:77 % 坐标差*4+1

for m=1:65

if d2((19*12+2)-(3*(p-1)+1),3*(m-1)+1)<0

high(p,m)=nan;

end

end

end

[X,Y]=meshgrid(105.75:.25:122,4:.25:23);

[m1 m2]=contourf(X,Y,high);

caxis([-30,30]);

colorbar

hold on

u = nan(77,66);

v = nan(77,66);

for i=2:76

k=29.75-0.25*(i-2);

f=2*(7.27e-5)*sin(k*(pi/180));

for j=2:65

g=1000;

v(i,j)=(1/f)*g*[high(i,j+1)-high(i,j-1)]/[(111*1/2)*cos(k*(pi/180))*(10^5)];

u(i,j)=-(1/f)*g*[high(i+1,j)-high(i-1,j)]/[(111*1/2)*(10^5)];

z(i,j)=[(v(i-1,j-1)^2)+(u(i-1,j-1)^2)]^(1/2);

end

xlabel('Lon[^oE]');

ylabel('Lat[^oN]');

end

legend('2cm/s',4);

set(gcf,'unit','centimeters','position',[10 2 19 16]);

quiver(X,Y,u,v,3)

hold off