Relatively simple problem with matrix and loops

2 views (last 30 days)
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

Accepted Answer

Fangjun Jiang
Fangjun Jiang on 14 Oct 2011
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
  4 Comments
Jules Ray
Jules Ray on 17 Oct 2011
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')

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!