clear;clc;
tic
[DEM,y]=geotiffread('截取dem.tif');
TREE=geotiffread('截取针叶林.tif');%12.tif是林地
info = geotiffinfo('截取dem.tif');
[M N]=size(TREE);
for j = 1:N
I1(1,j)=0;
I1(7,j)=0;
end
for i = 1:M
I1(i,1)=0;
I1(i,10)=0;
end
for i=2:M-1
for j=2:N-1
if TREE(i,j)==5
a1=(DEM(i-1,j-1)-DEM(i,j))/1.4;
a2=DEM(i,j-1)-DEM(i,j);
a3=DEM(i-1,j)-DEM(i,j);
a4=DEM(i+1,j+1)-DEM(i,j)/1.4;
a5=DEM(i+1,j-1)-DEM(i,j)/1.4;
a6=DEM(i-1,j+1)-DEM(i,j)/1.4;
a7=DEM(i+1,j)-DEM(i,j);
a8=DEM(i,j+1)-DEM(i,j);
if a8>=a7&&a8>=a6&&a8>=a5&&a8>=a4&&a8>=a3&&a8>=a2&&a8>=a1&&a8>=0&&TREE(i,j+1)==5
I1(i,j)=0;
elseif a7>=a8&&a7>=a6&&a7>=a5&&a7>=a4&&a7>=a3&&a7>=a2&&a7>=a1&&a7>=0&&TREE(i+1,j)==5
I1(i,j)=0;
elseif a6>=a8&&a6>=a7&&a6>=a5&&a6>=a4&&a6>=a3&&a6>=a2&&a6>=a1&&a6>=0&&TREE(i-1,j+1)==5
I1(i,j)=0;
elseif a5>=a8&&a5>=a7&&a5>=a6&&a5>=a4&&a5>=a3&&a5>=a2&&a5>=a1&&a5>=0&&TREE(i+1,j-1)==5
I1(i,j)=0;
elseif a3>=a8&&a3>=a7&&a3>=a5&&a3>=a4&&a3>=a6&&a3>=a2&&a3>=a1&&a3>=0&&TREE(i-1,j)==5
I1(i,j)=0;
elseif a2>=a8&&a2>=a7&&a2>=a5&&a2>=a4&&a2>=a3&&a2>=a6&&a2>=a1&&a2>=0&&TREE(i,j-1)==5
I1(i,j)=0;
elseif a1>=a8&&a1>=a7&&a1>=a5&&a1>=a4&&a1>=a3&&a1>=a2&&a1>=a6&&a1>=0&&TREE(i-1,j-1)==5
I1(i,j)=0;
else
I1(i,j)=1;
end
else
I1(i,j)=0;
end
end
end
geotiffwrite('result.tif',I1,y,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
subplot(1,3,1),imshow(DEM);
subplot(1,3,2),imshow(TREE);
subplot(1,3,3),imshow(I1);
toc
tic
[DEM,y]=geotiffread('截取dem.tif');
TREE=geotiffread('截取针叶林.tif');%12.tif是林地
info = geotiffinfo('截取dem.tif');
[M N]=size(TREE);
for j = 1:N
I1(1,j)=0;
I1(7,j)=0;
end
for i = 1:M
I1(i,1)=0;
I1(i,10)=0;
end
for i=2:M-1
for j=2:N-1
if TREE(i,j)==5
a1=(DEM(i-1,j-1)-DEM(i,j))/1.4;
a2=DEM(i,j-1)-DEM(i,j);
a3=DEM(i-1,j)-DEM(i,j);
a4=DEM(i+1,j+1)-DEM(i,j)/1.4;
a5=DEM(i+1,j-1)-DEM(i,j)/1.4;
a6=DEM(i-1,j+1)-DEM(i,j)/1.4;
a7=DEM(i+1,j)-DEM(i,j);
a8=DEM(i,j+1)-DEM(i,j);
if a8>=a7&&a8>=a6&&a8>=a5&&a8>=a4&&a8>=a3&&a8>=a2&&a8>=a1&&a8>=0&&TREE(i,j+1)==5
I1(i,j)=0;
elseif a7>=a8&&a7>=a6&&a7>=a5&&a7>=a4&&a7>=a3&&a7>=a2&&a7>=a1&&a7>=0&&TREE(i+1,j)==5
I1(i,j)=0;
elseif a6>=a8&&a6>=a7&&a6>=a5&&a6>=a4&&a6>=a3&&a6>=a2&&a6>=a1&&a6>=0&&TREE(i-1,j+1)==5
I1(i,j)=0;
elseif a5>=a8&&a5>=a7&&a5>=a6&&a5>=a4&&a5>=a3&&a5>=a2&&a5>=a1&&a5>=0&&TREE(i+1,j-1)==5
I1(i,j)=0;
elseif a3>=a8&&a3>=a7&&a3>=a5&&a3>=a4&&a3>=a6&&a3>=a2&&a3>=a1&&a3>=0&&TREE(i-1,j)==5
I1(i,j)=0;
elseif a2>=a8&&a2>=a7&&a2>=a5&&a2>=a4&&a2>=a3&&a2>=a6&&a2>=a1&&a2>=0&&TREE(i,j-1)==5
I1(i,j)=0;
elseif a1>=a8&&a1>=a7&&a1>=a5&&a1>=a4&&a1>=a3&&a1>=a2&&a1>=a6&&a1>=0&&TREE(i-1,j-1)==5
I1(i,j)=0;
else
I1(i,j)=1;
end
else
I1(i,j)=0;
end
end
end
geotiffwrite('result.tif',I1,y,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
subplot(1,3,1),imshow(DEM);
subplot(1,3,2),imshow(TREE);
subplot(1,3,3),imshow(I1);
toc