Cross sectional Profile of lakes by satellite image and polygon files?

2 views (last 30 days)
Hi there, I have a Shape file (about 1000 polygons that represent lakes) and a raster file (satellite image, geotiff) from the same area. I want to make cross sectional profiles (in the same direction, see the image attached) through each lake from the values of the raster (pixel values that the line crosses in each lake) and finally create a matrix which show each row as a lake and the numbers in each row that represents the value of cross section line on the raster. I know it won't be a matrix because the lakes have different sizes and the their cross sectional profile values are different. But any idea of how to do the process and finally store it in such a way is really appreciated. Thanks in advance
  2 Comments
John BG
John BG on 22 May 2018
Hi Marmar
is the angle for the cross sections constant, the one in the attached image, or do you need the angle to be a parameter?

Sign in to comment.

Answers (2)

Image Analyst
Image Analyst on 22 May 2018
To make profiles of image values along a line across the image, use improfile().
  1 Comment
Image Analyst
Image Analyst on 26 May 2018
Presumably he already knows where the lakes are. If he doesn't, and the images are really those computer graphics images, then you can simply do color segmentation to find the blue things. You can get the orientation, majorAxisLength and centroid from regionprops . Or use the Feret diameter from Steve's blog https://blogs.mathworks.com/steve/2018/04/17/feret-properties-wrapping-up/ or from my attached demo.
So once you have that you can easily find the endpoints of a line that goes through the blob, and then you can use improfile().
For real images, like from a satellite, you'd have to hyperspectral image segmentation to find the water and lakes. Probably not very hard.
I have no idea what he meant when he said this: "I want to apply the wind direction."

Sign in to comment.


N/A
N/A on 26 May 2018
Edited: John Kelly on 1 Jun 2018
1.-
For a given input angle, constant as you mentioned throughout the entire image, it's far easier rotating the image and then applying vertical rasters only.
Some time ago I wrote a function to rotate images given an angle and a pivot point; imrotatex.m attached to this answer, so the solution gets significantly simplified if instead of applying parallel raster lines of the input angle, we rotate the image the desired angle, and then apply, let's say, just vertical raster lines seeking the water.
The attached example script, explained
2.-
Assume the following
% convention 1: X axis: top horizontal side of image
% 2: Y axis: left vertical side of image
% 3: a [degree]: between bottom horizontal side of image and
% chosen cross-section line.
% 4: a range [0 180]
% 5: limit a resolution to da decimals da range [0 3]
% 6: sweeping always CCW
3.- acquiring the image
clear all;close all;clc
A=imread('001.png');
figure(1);imshow(A);
[sz1 sz2 sz3]=size(A)
4.- measuring size of image
x1=floor(sz2/2);y1=floor(sz1/2);
5.- Setting the raster angle
dlg1='Angle'; % choose angle
prompt1={'angle in degree [º]:'};numlines1=1;defaultans1={'30'};
angle1=str2num(cell2mat(inputdlg(prompt1,dlg1,numlines1,defaultans1)));
.
6.-
Choose -45º to align the single raster printed on the input image with the vertical axis.
Another convention:
negative angles; Counter Clock Wise
positive angles; Clock Wise
7.-
Rotating and showing the expanded image to be used
[A2,x_1,y_1,x_c,y_c,hf]=imrotatex(A,pi/2-angle1,x1,y1);
.
Ignore the red circle and the star, it's all data that shows only as reference, but in the next image it's omitted.
8.-
Now this is the start point
.
.
close all;
figure(1);imshow(A2)
[sz1_A2,sz2_A2,sz3_A2]=size(A2)
9.-
Define what is water with 3 colour ranges. The pixels that have all 3 RGB channels within these ranges are considered water.
range_R_blue=[45 75] % define water
range_G_blue=[90 135]
range_B_blue=[140 255]
10.-
main loop
L={}
for k=1:1:sz2_A2
L01=A2(:,k,1)
L02=A2(:,k,2)
L03=A2(:,k,3)
% hold on; plot(k*ones(1,sz2_A2),[1:1:sz2_A2],'r*') % test
% A2(:,k,1)=10;A2(:,k,2)=10;A2(:,k,3)=250 % another check
% figure(1);imshow(A2)
nL01_1=find(L01>=range_R_blue(1))
nL01_2=find(L01<=range_R_blue(2))
nL01=intersect(nL01_1,nL01_2)
nL02_1=find(L02>=range_G_blue(1))
nL02_2=find(L02<=range_G_blue(2))
nL02=intersect(nL02_1,nL02_2)
nL03_1=find(L03>=range_B_blue(1))
nL03_2=find(L03<=range_B_blue(2))
nL03=intersect(nL03_1,nL03_2)
nL=intersect(nL01,intersect(nL02,nL03))
hold on
plot(k*ones(1,numel(nL)),nL,'go')
% stem(nL)
L{k}=nL
% 003
end
save lake_search_vars.mat
11.-
Note the commented lines
% hold on; plot(k*ones(1,sz2_A2),[1:1:sz2_A2],'r*') % test
% A2(:,k,1)=10;A2(:,k,2)=10;A2(:,k,3)=250 % another check
% figure(1);imshow(A2)
are useful to check whether a single raster is detecting the right pixels. For instance for a given
k=256
hold on; plot(k*ones(1,sz2_A2),[1:1:sz2_A2],'r*') % test
A2(:,k,1)=10;A2(:,k,2)=10;A2(:,k,3)=250 % another check
figure(1);imshow(A2)
the sought matrix containing the coordinates of water is the cell variable
L
contained in the saved workspace
lake_search_vars.mat
also attached to this answer
load .mat file with
load lake_search_vars.mat
Address a single raster segment
k=20;L{k}
.
.
Additional comments
.
1.-
The choice of what's water and what's not, through the 3 ranges above mentioned, is important.
I Just checked a few values, so I guessed that for a pixel to be considered blue water it had to have the values chosen above.
With such ranges, the result is:
2.-
This imrotatex.m has slight changes compared to the original imrotatex.m that solved the following questions:
.
.
3.-
Following, copy of the attached support function
function imrotatex.m
function [B4,x1,y1,xc,yc,ax]=imrotatex(A,angle,x1,y1)
% ithis function rotates image A around Z axis centered on [x1 y1]
% returning
% B4: rotated image .
% [x1 y1] on original [0 0] reference, any point within image A.
% [xc yc] centre point of original image, according to initial [0 0].
%
% in version 1 point [x1 y1] is introduced with mouse.
%
% angle defined in degree [0 360] if >0 then CCW if CCW< rotation CW.
% angle fractions in decimal: 360.50 would be 360º30'0''.
%
% imrotatex puts:
% original image centre with o
% mouse marked rotation point with *
% circle around image to show
%
% depending upon in what sector [x1 y1] is defined, padarray needs different options
% in the following sketch if o is the [xc yc] center of the original image.
% The sectors are defined as follows
% [0 0] ------------> x
% | 1 2
% | o
% | 4 3
% y
%
% case 1: x1<xc & y1<yc
% case 2: x1>cx & y1<yx
% case 3: x1>xc & y1>yc
% case 4: x1<xc & y1>yc
%
% The resulting image is saved to .jpeg file
%
% if attempting to crop, imtransform warns that the resulting image in the
% result figure has been soomed out to 67%. If you intend to crop the
% rotated image, you may want to first zoom in to have same scale factor
% as original figure
%
% included
% test file: example_imrotatex.m
% sample image: countryside.jpg
% only tested with .jpg images
%%%%%%%%%%%%%%%%%%
% version 1. June 3rd 2016
% Author: John Bofarull Guix,
% feedbacl. comments and suggestions: jgb2012@sky.com
%
%%%%%%%%%%%%%%%%%%
N=255; % value used in padarray fill-up
if abs(angle)<1
disp('choose angle > 1º');B4=A;
return; end;
[sy sx sz]=size(A); % measure image size
xc=floor(sx/2);yc=floor(sy/2); % find image centre [xc yc]
hold on;figure(1);plot(x1,y1,'r*','MarkerSize',10); % visualize both [xc yc] and [x1 y1]
hold on;figure(1);plot(xc,yc,'ro','MarkerSize',10);
if (x1==xc) % make sure x1<>xc and y1<>yc
x1=x1+2*(randi([0 1],1,1)-.5); end;
if (y1==yc)
y1=y1+2*(randi([0 1],1,1)-.5); end;
ks=0;
if (x1<xc & y1<yc) ks=1; end;
if (x1>xc & y1<yc) ks=2; end;
if (x1>xc & y1>yc) ks=3; end;
if (x1<xc & y1>yc) ks=4; end;
t=0:2*pi/100:2*pi;
switch ks
case 1 % x1<xc & y1<yc
D=norm([x1 y1]-[sx sy]);
dx=D+x1-sx; % 1st pair of padding margins
alpha=atand((sx-x1)/(sy-y1));
s=(2*D^2-2*D^2*cosd(alpha))^.5;
dy=(s^2-(sx-x1)^2)^.5;
% dy=D+y1
dx=floor(dx);dy=floor(dy);
B=padarray(A,[dy,dx],N,'both');
figure(2);imshow(B);
[sby sbx sbz]=size(B);
hold on;figure(2);plot(x1+dx,y1+dy,'r*','MarkerSize',10);
hold on;figure(2);plot(xc+dx,yc+dy,'ro','MarkerSize',10);
% the origin has been displaced [-dx -dy]. the point [x1 y1] to become rotation Z axis now has coordinates [dx+x1 +dy+y1]
% visualizing the outer rotation circle
x0d=D*cos(t); xd=x0d+dx+x1;
y0d=D*sin(t); yd=y0d+dy+y1;
hold all;figure(2);plot(xd,yd,'r','LineWidth',1.5);
% visualizing axes after 2nd padarray and origin displacement
hold all;figure(2);plot([1:1:sbx],ones(1,sbx).*(dy+y1),'w','LineWidth',1.5); % new x axis
hold all;figure(2);plot(ones(1,sby).*(dx+x1),[1:1:sby],'w','LineWidth',1.5); % new y axis
dx2=2*D-sbx;dx2=floor(dx2); % calculating 2nd padarray margins
dy2=2*D-sby;dy2=floor(dy2);
B3=padarray(B,[dy2 dx2],N,'pre');
[sb3y sb3x sb3z]=size(B3);
figure(4);imshow(B3);
hold all;figure(4);plot(x1+dx+dx2,y1+dy+dy2,'r*','MarkerSize',10); % visualizing initial image centre [xc yc] and and chosen rotation point [x1 y1]
hold all;figure(4);plot(xc+dx+dx2,yc+dy+dy2,'ro','MarkerSize',10);
xd2=xd+dx2; % visualizing outer circle after 2nd padarray
yd2=yd+dy2;
hold all;figure(4);plot(xd2,yd2,'r','LineWidth',2.5);
% visualizing axes after 2nd padarray and origin displacement
hold all;figure(4);plot([1:1:sb3x],ones(1,sb3x).*(dy+y1+dy2),'w','LineWidth',1.5); % new x axis
hold all;figure(4);plot(ones(1,sb3y).*(dx+x1+dx2),[1:1:sb3y],'w','LineWidth',1.5); % new y axis
ax=gca;
case 2 % x1>cx & y1<yx
D=norm([x1 y1]-[0 sy]);
dx=D-x1; % 1st pair of padding margins
alpha=atand((sy-y1)/(sx-x1));
s=(2*D^2-2*D^2*cosd(alpha))^.5;
dy=(s^2-(sx-x1)^2)^.5;
dy=D-y1;
dx=floor(dx);dy=floor(dy);
B=padarray(A,[dy,dx],N,'pre');
figure(2);imshow(B);
[sby sbx sbz]=size(B);
hold on;figure(2);plot(x1+dx,y1+dy,'r*','MarkerSize',10);
hold on;figure(2);plot(xc+dx,yc+dy,'ro','MarkerSize',10);
% the origin has been displaced [-dx -dy]. the point [x1 y1] to become rotation Z axis now has coordinates [dx+x1 +dy+y1]
% visualizing the outer rotation circle
x0d=D*cos(t); xd=x0d+dx+x1;y0d=D*sin(t); yd=y0d+dy+y1;
hold all;figure(2);plot(xd,yd,'r','LineWidth',1.5);
% visualizing axes after 2nd padarray and origin displacement
hold all;figure(2);plot([1:1:sbx],ones(1,sbx).*(dy+y1),'w','LineWidth',1.5); % new x axis
hold all;figure(2);plot(ones(1,sby).*(dx+x1),[1:1:sby],'w','LineWidth',1.5); % new y axis
dx2=D-(sx-x1);dx2=floor(dx2); % calculating 2nd padarray margins
dy2=D-(sy-y1);dy2=floor(dy2);
B3=padarray(B,[dy2 dx2],N,'post');
[sb3y sb3x sb3z]=size(B3);
figure(4);imshow(B3);
xc2=floor(sb3x/2);yc2=floor(sb3y/2);
hold all;figure(4);plot(xc2,yc2,'r*','MarkerSize',10); % visualizing initial image centre [xc yc] and and chosen rotation point [x1 y1]
hold all;figure(4);plot(xc2-(x1-xc),yc2-(y1-yc),'ro','MarkerSize',10);
xd2=xd+xc2; yd2=yd+yc2; % visualizing outer circle after 2nd padarray
hold all;figure(4);plot(xd,yd,'r','LineWidth',2.5);
% visualizing axes after 2nd padarray and origin displacement
hold all;figure(4);plot([1:1:sb3x],ones(1,sb3x).*yc2,'w','LineWidth',1.5); % new x axis
hold all;figure(4);plot(ones(1,sb3y).*xc2,[1:1:sb3y],'w','LineWidth',1.5); % new y axis
ax=gca;
case 3 % x1>xc & y1>yc
D=norm([x1 y1]-[0 0]);
dx=D-x1;dy=D-y1; % 1st pair of padding margins
dx=floor(dx);dy=floor(dy);
B=padarray(A,[dy dx],N,'pre');
figure(2);imshow(B);
[sby sbx sbz]=size(B);
hold on;figure(2);plot(x1+dx,y1+dy,'r*','MarkerSize',10);
hold on;figure(2);plot(xc+dx,yc+dy,'ro','MarkerSize',10);
x0d=D*cos(t); xd=x0d+dx+x1;y0d=D*sin(t); yd=y0d+dy+y1;
hold all;figure(2);plot(xd,yd,'r','LineWidth',1.5);
% visualizing axes after 2nd padarray and origin displacement
hold all;figure(2);plot([1:1:sbx],ones(1,sbx).*(dy+y1),'w','LineWidth',1.5); % new x axis
hold all;figure(2);plot(ones(1,sby).*(dx+x1),[1:1:sby],'w','LineWidth',1.5); % new y axis
dx2=2*D-sbx;dx2=floor(dx2); % calculating 2nd padarray margins
dy2=2*D-sby;dy2=floor(dy2);
B2=padarray(B,[0 dx2],N,'post');
figure(3);imshow(B2);
B3=padarray(B2,[dy2 0],N,'post');
[sb3y sb3x sb3z]=size(B3)
figure(4);imshow(B3);
xc2=floor(sb3x/2);yc2=floor(sb3y/2);
% visualizing initial image centre [xc yc] and and chosen rotation point [x1 y1]
hold all;figure(4);plot(xc2,yc2,'r*','MarkerSize',10)'
hold all;figure(4);plot(xc2-(x1-xc),yc2-(y1-yc),'ro','MarkerSize',10)'
xd2=x0d+xc2;yd2=y0d+yc2'; % visualizing outer circle after 2nd padarray
hold all;figure(4);plot(xd2,yd2,'r','LineWidth',2.5)'
% visualizing axes after 2nd padarray and origin displacement
hold all;figure(4);plot([1:1:sb3x],ones(1,sb3x).*yc2,'w','LineWidth',1.5); % new x axis
hold all;figure(4);plot(ones(1,sb3y)*xc2,[1:1:sb3y],'w','LineWidth',1.5); % new y axis
ax=gca;
case 4 % x1<xc & y1>yc
D=norm([x1 y1]-[sx 0]); % D is distance between [x1 y1] and [xc yc], also D=pdist2([x1 y1],[sx 0])
dx=D+x1-sx; % 1st pair of padding margins
alpha=atand((sx-x1)/y1)
s=(2*D^2-2*D^2*cosd(alpha))^.5
dy=(s^2-(sx-x1)^2)^.5
dx=floor(dx);dy=floor(dy);
B=padarray(A,[dy,dx],N,'both');
figure(2);imshow(B)
[sby sbx sbz]=size(B)
hold on;figure(2);plot(x1+dx,y1+dy,'r*','MarkerSize',10)
hold on;figure(2);plot(xc+dx,yc+dy,'ro','MarkerSize',10)
% the origin has been displaced [-dx -dy]. the point [x1 y1] to become rotation Z axis now has coordinates [dx+x1 +dy+y1]
% visualizing the outer rotation circle
x0d=D*cos(t); xd=x0d+dx+x1;y0d=D*sin(t); yd=y0d+dy+y1;
hold all;figure(2);plot(xd,yd,'r','LineWidth',1.5);
% visualizing axes after 2nd padarray and origin displacement
hold all;figure(2);plot([1:1:sbx],ones(1,sbx).*(dy+y1),'w','LineWidth',1.5) % new x axis
hold all;figure(2);plot(ones(1,sby).*(dx+x1),[1:1:sby],'w','LineWidth',1.5) % new y axis
dx2=2*D-sbx;dx2=floor(dx2); % calculating 2nd padarray margins
dy2=2*D-sby;dy2=floor(dy2);
B2=padarray(B,[0 dx2],N,'pre');
figure(3);imshow(B2);
B3=padarray(B2,[dy2 0],N,'post');
[sb3y sb3x sb3z]=size(B3)
figure(4);imshow(B3);
% visualizing initial image centre [xc yc] and and chosen rotation point [x1 y1]
hold all;figure(4);plot(x1+dx+dx2,y1+dy,'r*','MarkerSize',10)
hold all;figure(4);plot(xc+dx+dx2,yc+dy,'ro','MarkerSize',10)
% visualizing outer circle after 2nd padarray
xd2=xd+dx2;yd2=yd;
hold all;figure(4);plot(xd2,yd2,'r','LineWidth',2.5)
% visualizing axes after 2nd padarray and origin displacement
hold all;figure(4);plot([1:1:sb3x],ones(1,sb3x).*(dy+y1),'w','LineWidth',1.5) % new x axis
hold all;figure(4);plot(ones(1,sb3y).*(dx+x1+dx2),[1:1:sb3y],'w','LineWidth',1.5) % new y axis
ax=gca;
otherwise
disp('\nerror\n');
return
end
%
% imsave(B3)
imwrite(B3,'padded_image.jpg');
% padded image: B3
% size square containing padded image: [sb3x sby3]
% center padded image: [xc2 yc2]
% points outer circle centre before rotation [xd2 yd2]
% radius: D
B4=imrotate(B3,angle);
[sb4y sb4x sb4z]=size(B4);
figure(5);imshow(B4);
% visualizing initial image centre [xc yc] and and chosen rotation point [x1 y1]
hold all;figure(5);plot(floor(sb4x/2),floor(sb4y/2),'r*','MarkerSize',15);
% hold all;figure(5);plot(xc+dx+dx2,yc+dy,'ro','MarkerSize',10);
% visualizing outer circle after 2nd padarray
xd4=x0d+floor(sb4x/2);yd4=y0d+floor(sb4y/2);
hold all;figure(5);plot(xd4,yd4,'r','LineWidth',2.5);
if angle>0 str1='CCW'; end;
if angle<0 str1='CW'; end;
figure(5);legend(['angle : ' num2str(angle) ' º ' str1],['Radius : ' num2str(floor(D)) ' pixels'],'Location','Northwest');
imwrite(B4,'im_padded_rotated.jpg');
% clear alpha dx dx2 dy2 ks sbx sby sbz sx sy sz t xc xd xd2 yc yc2 yd; % clean as you go
clc;
end

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!