Relatively simple problem with matrix and loops

 i got a problem with a script and is described next:

this is the script that i'm preparing. Is focused to clip a topography in several parts (Topoclip1, Topoclip2... Topoclipi) based on polygons included in a shapefile.

The main problem is the loop, i dont know what's happen when i try to run the scripts the results is (??? Subscripted assignment dimension mismatch.)

 %%Inputs

 clear
 addpath /Users/mac/Documents/MATLAB/intersection/stat1
 stationsdir='/Users/mac/Documents/MATLAB/intersection/stat1';

 a=load('topo.txt');  %load the raster (previously grided xyz table)

 cd(stationsdir)
 %cd(station)
 station='stat'; 
 file=sprintf('%s_clip.shp',station); 
 p=shaperead(file);  
 nim=numel(p);

 % Extract rectangular profile corners for each profile

 for i=1:nim
    data=p(i,1); %for the current extention of data...
    x=data.X; %extract x values from shapefile
    x=x'; %transpose elements
    x(end)=[];
    y=data.Y; %extract y values from shapefile
    y=y'; %transpose elements
    y(end)=[]; 
    xy=[x y]; %conform the resulting matrix
    xy(end,:)=[]; %delete the redundating last row and fix a 4 vertices polygon

   %to here all works ok... the problems is placed somewhere in the next lines

    taga=num2str(i);
    v = genvarname(['Prof' taga]);
    eval([v '=xy']);

    Xp=xy(:,1); %rectangular profile X coordinate
    Yp=xy(:,2); %rectangular profile Y coordinate
    Xt= a(:,1); %Lon1
    Yt= a(:,2); %Lat1
    h= a(:,3);

    %in = inpolygon(Xt,Yt,Xp,Yp);
    %plot(Xp,Yp,x(in),Yt(in),'r+')

    [loni,lati] = polybool('intersection',Xt,Yt,Xp,Yp);
    [lati loni];

    A(:,1)=a(:,1);
    A(:,2)=a(:,2);
    A(:,3)=a(:,3);
    B(:,1)=loni';
    B(:,2)=lati';

    out=A(ismember(A(:,1:2),B,'rows'),:);

    taga=num2str(i);
    v = genvarname(['Topoclip' taga]);
    eval([v '=out']);

    %scatter(Topoclip1(:,1), Topoclip1(:,2),20,Topoclip1(:,3),'filled')

end

the ides is to obtain an i number of Topoclip (ie: Topoclip1, Topoclip2, Topocli3....) eachone formed by the clipped topography.

Please give me a hand.... if you can...

thanx

Answer 1

[Fangjun Jiang]

If you run the following, you have the same problem. The error usually tells you which line of the code did the error happen. So check the assignment inside the loop to see how it happens. You can put a break point in the code and run the code line by line

a=rand(3);
a(:,1)=1:5

[Jules Ray]

Thanks Fangjun, i think i have found something...

Specifically about this function:

in = inpolygon(Xt,Yt,Xp,Yp);

te result (in) is a logical value, so, how do you write a code to obtain the X,Y values (xy coordinates) of the points inside the polygon.... i think after fix this trouble the loop can run good

[Fangjun Jiang]

Sorry, I am not familiar with polygon. If you could make your code runnable for others without your path and data, I might be able to find where is the error. In fact, if you follow the error message or add break point in your code, you could easily find out.

[Walter Roberson]

Jules,

After your line

in = inpolygon(Xt,Yt,Xp,Yp);

then the points that are inside the polygon are Xt(in) and Yt(in)

If you need the indices of these for some purpose, then they are

find(in)

[Jules Ray]

Thanx Walter, your solution has let me finish the code, and now works perfectly. I would like to share this with everyone as a demonstration of gratitude:

%INTER sript to multiclip an geotiff topography, using .shp polygons

%By Julius J.

%a .xyz topography

%a shapefile containing rectangular polygons (rectangular profiles)

%% Inputs

clear

addpath /Users/mac/Documents/MATLAB/intersection/stat1

stationsdir='/Users/mac/Documents/MATLAB/intersection/stat1';

a=load('topo.txt'); %load the raster

%a=geotiffread('topotif.tif') %in caseo of a geotiff input

cd(stationsdir)

%cd(station)

station='stat';

file=sprintf('%s_clip.shp',station); %load the .shp rectangular profiles

p=shaperead(file);

nim=numel(p);

%% Extract corners for each rectangular profile and clip with topography

for i=1:nim

data=p(i,1); %for the current extention of data...

x=data.X; %extract x values from shapefile

x=(x'); %transpose elements

x(end)=[];

y=data.Y; %extract y values from shapefile

y=(y'); %transpose elements

y(end)=[];

xy=[x y]; %conform the resulting matrix

xy(end,:)=[]; %delete the redundating last row and fix a 4 vertices polygon

taga=num2str(i); %assign names and vertice coordinates for each polygon

v = genvarname(['Prof' taga]);

eval([v '=xy']);

Xp=xy(:,1); %rectangular profile X lon2

Yp=xy(:,2); %rectangular profile Y lat2

Xt= a(:,1); %Lon1

Yt= a(:,2); %Lat1

h= a(:,3); %still not used z value of topography

in = inpolygon(Xt,Yt,Xp,Yp); %clip only X Y data of topography

loni=Xt(in);

lati=Yt(in); ex=[loni lati]; %ex is defined as the output of the cliped xy data

ex(end,:)=[];

A(:,1)=a(:,1); %we redefine A as the topography data

A(:,2)=a(:,2);

A(:,3)=a(:,3);

out=A(ismember(A(:,1:2),ex,'rows'),:); %add the z value to clipped data, generating the final table

taga=num2str(i);

v = genvarname(['Topoclip' taga]); %Asign names for each cliped area

eval([v '=out']); %Asign xyz values for each cliped area

end

clear data Xp Xt Yp Yt h lati loni xy x y p taga a nim out ex i

%% Graphic outputs

hold on

scatter(A(:,1), A(:,2),A(:,3),'o')

scatter(Topoclip1(:,1), Topoclip1(:,2),20,Topoclip1(:,3),'filled')

scatter(Topoclip2(:,1), Topoclip2(:,2),20,Topoclip2(:,3),'filled')

scatter(Topoclip3(:,1), Topoclip3(:,2),20,Topoclip3(:,3),'filled')