Matlab 根據 shp 裁剪矩陣/圖像 函數

1、全代碼

function varargout=tailorsheng(varargin)
%% 此函數用於根據 shp 裁剪矩陣/圖像
% 輸入：
%       P2file shpfile
%       TDMSPfile tiffile，可以不是tif
%       str  要裁剪誰
%       mode 1省0國家
% 輸出：
%       sheng shp
%       photo mat圖像
%       ex    經緯度極值
% 調用：
%       TDMSPfile='D:\復習資料\自學\7_科研\夜光遙感\data\DMSP\F121999.v4\F121999.v4b_web.stable_lights.avg_vis.tif';
%       P2file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4p.shp';%省界多邊形
%       str='北京市';
%       mode=1;
%       [sheng,photo,ex]=tailorsheng(P2file,TDMSPfile,str);
%       [sheng,photo,ex]=tailorsheng(P2file,TDMSPfile,str,mode);
%-------------------------------------------------------------------
    %%%%    Authors:   Bill O'Hanlon
    %%%%    EMAIL:     [e-mail]
    %%%%    DATE:      24-08-2020
%% 輸入
disp('-------function tailorsheng------');
tic %開始計時
mode=1;
if nargin==3
    P2file=varargin{1}; %注意，這里是{, (會是元胞
    TDMSPfile=varargin{2};
    str=varargin{3};
elseif nargin==4
    P2file=varargin{1}; %注意，這里是{, (會是元胞
    TDMSPfile=varargin{2};
    str=varargin{3};
    mode=varargin{4};%這里不做變換的話會是元胞
else
    disp('參數不合要求！');
    return;
end
%% 第一步，提取目標省份的 shp，注意看其范圍
[henan,ex]=drawsheng(P2file,str,mode);
disp('shp文件搞定！');
%% 第二步，把圖像大致范圍剪出來，這步根據圖像種類不同，代碼也不一樣，Attention!!
TDMSP=imread(TDMSPfile); %DMSP范圍 180°W-180°E，65°S-75°N
Xmax=ex(1);              %X是經度，Y是緯度
Ymax=ex(2);              %這些已經取整了。
Xmin=ex(3);
Ymin=ex(4);
Xmax1=Xmax+180;          %這是裁剪時用的
Ymax1=75-Ymax;
Xmin1=Xmin+180;
Ymin1=75-Ymin;
Henan=TDMSP(Ymax1*120:Ymin1*120,Xmin1*120:Xmax1*120);%根據自己需要裁減
disp('粗裁剪完成！');
%% 第三步，計算不規則邊界的內切圓和外接圓
Y=Ymin:0.0083333333:Ymax; %分辨率0.0083333333°  1°=120
X=Xmin:0.0083333333:Xmax; %求邏輯矩陣用到
[a,b]=size(Henan);
Y2=repmat(Y',1,b);
X2=repmat(X,a,1);
Y2=flipud(Y2);                     %這個緯度得上下翻轉一下。
xv=henan.X;yv=henan.Y;             %提取邊界
xv = xv(1:end-1); yv = yv(1:end-1);%把最后一個NaN去掉
disp('需要計算邏輯的區域為 Figure 2 紅色區域');
bianjie=[xv;yv];
[zhongxin1,zhongxin2,smallR,bigR]=getZhongxin(bianjie,X2,Y2);
%% 第四步，根據內切圓和外接圓+邊界進行裁剪
disp('開始計算邏輯矩陣..');
photo=zeros(a,b);
photo(:)=nan;
for i=1:a
    for j=1:b
        if mod(i*b+j,10000)==0
            disp([int2str(i*b+j),' / ',int2str(a*b), ' OK!']);
        end
        x=X2(i,j);
        y=Y2(i,j);
        if norm(zhongxin1-[x;y])<=smallR
            photo(i,j)=Henan(i,j);%內接圓里面的都在
            continue;
        elseif norm(zhongxin2-[x;y])>=bigR
            continue;%外接圓外面的都不在
        elseif inpolygon(x,y,xv,yv)
            photo(i,j)=Henan(i,j);
        end
    end
end
disp('細裁剪完成！');
%% 第五步，畫圖，裁剪后的
[x,x1,y,y1] = getxy(X,Y);
x=(x-Xmin).*120; 
y=(y-Ymin).*120;
photo=flipud(photo);%contourf 上下翻轉一下，才變成imshow
figure;
contourf(photo,'LineStyle','none');
colormap(gray);colorbar %jet
set(gca,'XTick',x,'XTicklabel',x1);   %設置x,y軸
set(gca,'YTick',y,'YTicklabel',y1);
title([str,'夜光遙感影像']);
%% 輸出
if nargout==3
    varargout{1}=henan; 
    varargout{2}=photo;
    varargout{3}=ex;
elseif nargout==2
    varargout{1}=henan; 
    varargout{2}=photo;
elseif nargout==1
    varargout{1}=henan; 
elseif nargout==0
    return;
else
    disp('參數不合要求！');
    return;
end
disp('--------Finished!--------');
toc  %展示運行時間
end

依賴：

內切圓和外接圓 :https://blog.csdn.net/Gou_Hailong/article/details/108206335
扣shp：https://blog.csdn.net/Gou_Hailong/article/details/108209395
縮放矩陣或圖像：https://blog.csdn.net/Gou_Hailong/article/details/108206521
畫地圖注釋：https://blog.csdn.net/Gou_Hailong/article/details/108208442

注：這個函數裁剪小點的邊界還好，如果邊界太大或者圖像太大，運行時間會超級長的，這酸爽被我記錄在了下面博客中：

https://blog.csdn.net/Gou_Hailong/article/details/108147268

2、調用

調用代碼：

clc
clear
TDMSPfile='D:\復習資料\自學\7_科研\夜光遙感\data\DMSP\F121999.v4\F121999.v4b_web.stable_lights.avg_vis.tif';
P2file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4p.shp';%省界多邊形
str='北京市';
mode=1;
[sheng,photo,ex]=tailorsheng(P2file,TDMSPfile,str,mode);

結果：