|
关于这个函数，我最近（2020.11.1）又修改了一下，之前的版本第二步通用性不强，这波应该可以通用了。现将几个版本命个名： 2020-8-24 V1.0 有内切圆外切圆，第二步不通用 2020-11-1 V2.0 有内切圆外切圆，第二步通用 2020-11-1 V2.1 无内切圆外切圆，第二步通用
1、全代码
1.1 V1.0
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])=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
1.2 V2.0
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: 01-11-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范围 73.49781527°E-135.0894573°E，3.83544513°N-53.56042524°N
Xmax=ex(1); %X是经度，Y是纬度
Ymax=ex(2); %这些已经取整了。
Xmin=ex(3);
Ymin=ex(4);
Info=geotiffinfo(TDMSPfile); %读取图像信息，比如经纬度，分辨率
global PixelScale; %定义了几个全局变量，
global Pixel_;
global XBound;
global YBound;
PixelScale=Info.PixelScale(1,1);
Pixel_=1/PixelScale;
XBound=Info.BoundingBox(1,1);
YBound=Info.BoundingBox(2,2);
[a,b]=size(TDMSP);
if(Xmax>Info.BoundingBox(2,1)) Xmax=Info.BoundingBox(2,1); end
if(XminInfo.BoundingBox(2,2)) Ymax=Info.BoundingBox(2,2); end
if(YminInfo.BoundingBox(2,1)) Xmax=Info.BoundingBox(2,1); end
if(XminInfo.BoundingBox(2,2)) Ymax=Info.BoundingBox(2,2); end
if(Ymin |