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:     [email protected]
    %%%%    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);

结果：