- 积分
- 156
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2020-3-19
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
这是一个青藏高原地区降水图的脚本,现在想在画图时把shp文件用上,脚本如下(主要部分已用颜色突出),出图时体现不出shp文件信息,就大佬指点
%%
clear, clc
% matlabpool local 1
StartYear = 1979;
EndYear = 2014;
SpendYear = EndYear - StartYear + 1;
Level = 10;
DongTuJudge = zeros(34,18,SpendYear);
Deep = [0.0071,0.02792,0.06226,0.11887,0.21219,0.36607,0.61976,1.03803,1.72764,2.86461,4.73916,7.82977,12.92532,21.32647,35.17762];
ZeroDeep = zeros(34,18,12);
ZerosDeepYear = zeros(34,18,SpendYear);
% 调试变量
Surf1 = zeros(34,18,SpendYear);
Surf2 = zeros(34,18,SpendYear);
% 这部分是画图用的
Boundary = load('TP_boundary.txt');
VarName1 = Boundary(:,1);
VarName2 = Boundary(:,2);
for Year = StartYear:EndYear
ID = ['GLDAS_CLM_TP.',num2str(Year),'.nc'];
ORGANIC1= ncread(ID,'var85');
ORGANIC=ORGANIC1(:,:,:,:);
LAT = ncread(ID,'lat');
LON =ncread(ID,'lon');
SoldTemp = ORGANIC - 273.15*ones(size(ORGANIC));
year = Year-StartYear+1;
Winter_surf =(Winter(:,:,1)+Winter(:,:,3))/2;
Winter_shallow = (Winter(:,:,3)+Winter(:,:,4)+Winter(:,:,5))/3;
Winter_deep = (Winter(:,:,6)+Winter(:,:,8))/2;
MeanWinter_surf(:,:,year) = Winter_surf;
MeanWinter_shallow(:,:,year) = Winter_shallow;
MeanWinter_deep(:,:,year) = Winter_deep;
end
for i = 1:34
for j = 1:18
[i,j]
in = inpolygon(LON(i),LAT(j),VarName1,VarName2);
if in == 0%在边界范围外
MeanWinter_surf(i,j,:) = NaN;
MeanWinter_shallow(i,j,:) = NaN;
MeanWinter_deep(i,j,:) = NaN;
end
end
end
% 算趋势
dt=1;
z = zeros(34,18);
sl = zeros(34,18);
lcl = zeros(34,18);
ucl = zeros(34,18);
if SpendYear > 10
for i = 1:34
for j = 1:18
YY = zeros(SpendYear,1);
for k = 1:SpendYear
YY(k,1) = MeanWinter_shallow(i,j,k);
end
if ~isnan(max(YY)) && (max(YY) == min(YY) == 0)
[z(i,j), sl(i,j), lcl(i,j), ucl(i,j) ] =trendMMK( YY, dt );
end
if isnan(max(YY))
z(i,j) = NaN;
sl(i,j) = NaN;
lcl(i,j) = NaN;
ucl(i,j) = NaN;
end
end
end
else
disp('显著性检验')
end
minLON = min([min(VarName1),min(LON)]);
maxLON = max([max(VarName1),max(LON)]);
LonRange = [floor(minLON),ceil(maxLON)];
LonRange = double(LonRange);
minLAT = min([min(VarName2),min(LAT)]);
maxLAT = max([max(VarName2),max(LAT)]);
LatRange = [floor(minLAT),ceil(maxLAT)];
LatRange = double(LatRange);
LAT = double(LAT);
LON = double(LON);
shading flat;
axis equal;
axis([minLON maxLON minLAT maxLAT]);
cmap = colormap;
map=shaperead('2008_2010.shp');
geoshow(LAT,LON,[0 0 1],cmap);%把地图画出来
hold on;
hold on;
plot(VarName1,VarName2);
hold off;
contourf(LON,LAT,sl');
title('1960/10的风场')
|
|