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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 6947|回复: 2

[程序设计] 关于MATLAB用shp文件作图叠加的求助

[复制链接]

新浪微博达人勋

发表于 2020-12-5 12:32:09 | 显示全部楼层 |阅读模式

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

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

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的风场')






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

新浪微博达人勋

 成长值: 32430
发表于 2020-12-5 18:05:36 | 显示全部楼层

回帖奖励 +5 金钱

建议用m_map工具箱,m_shaperead和m_line
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-12-13 22:27:29 | 显示全部楼层
真的好巧
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

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

本版积分规则

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

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

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