How do I extract points from multiple circles?

I have plot circle according to the grid points and I have a problem to extract the points within each circle. This is my coding..

How should I extract the points from each circle. and save it in separate files according to the grid points. anyone can help?

clc;clear;

% Number of Points
numPoints = 5475;
data = load('Sample.txt');
x = data (:,2);
y = data (:,1);
z = data (:,3);%                      
plot(x, y, '.r');

% Make circle.
grid = load('Grid025.asc');
[row,column]=size(grid);
x0 = grid (:,2);
y0 = grid (:,1);
R = 0.09;


% Plot circle
for i=1:row
pos = [x0(i)-R, y0(i)-R, 2*R, 2*R];
rectangle('Position', pos, 'Curvature',[1 1]);
hold on;

% Determine how many points are in the circle .
count(i) = sum(((x-x0(i)).^2+(y-y0(i)).^2<=R^2));
listC = count.';

%Extract Data in Each Circle
distances = sqrt((x-x0(i)).^ 2 + (y-y0(i)).^ 2);
InsideCircle = distances <=R;
xS = x(InsideCircle); % Selected Lon
yS = y(InsideCircle); % Selected Lat
zS = z(InsideCircle); % Selected h
end

Answer 1

[Cedric]

We cannot test without your data, but the approach seems fine. You could use the PDIST2 function and extract everything in one shot, but as a first approach this is fine.

I will assume that the question is about saving the points, so here is small update of your code that should (if I understand your example well) generate an Excel file with 4 columns that store Circle ID, Point ID, Point X, Point Y, Point Z:

% Prealloc array of counts, and cell array for storing relevant point parameters per circle.
counts = zeros( row, 1 ) ;
output = cell( row, 1 ) ;

% Initialize figure.
figure() ;
hold on ;

% Process circles.
for circleId = 1 : row
    % Plot circle.
    pos = [x0(circleId)-R, y0(circleId)-R, 2*R, 2*R] ;
    rectangle( 'Position', pos, 'Curvature', [1 1] ) ;
    
    % Compute distance from all points to current circle center.
    distances = sqrt( (x - x0(circleId)).^2 + (y - y0(circleId)).^2 ) ;
    
    % Flag and count relevant points.
    isInsideCircle = distances <= R ;
    counts(circleId) = sum( isInsideCircle ) ;
    
    % Extract relevant parameters.
    pointsId = find( isInsideCircle ) ;
    pointsX = x(isInsideCircle) ; % Selected Lon
    pointsY = y(isInsideCircle) ; % Selected Lat
    pointsZ = z(isInsideCircle) ; % Selected h
    
    % Store in relevant output cell. We build the first column as a vector of elements equal
    % to circleId, whose size matches the size of the other columns (e.g. pointsId).
    output{circleId} = [circleId * ones( size( pointsId )), pointsId, pointsX, pointsY, pointsZ] ;
end

% Build data to export to Excel by concatenating all cells of the output cell array vertically
% into a large numeric array, converting it into a cell array, and concatenating to a first
% "header" row.
output = num2cell( vertcat( output{:} )) ;
output = [{'CircleID', 'PointID', 'PointX', 'PointY', 'PointZ'}; output] ;
 
% Export to Excel.
xlswrite( 'PointsPerCircle.xlsx', output ) ;

Note that I didn't test it, but it illustrates a possible approach. Let me know if you have any question.

[hanif hamden]

Your coding works perfectly but the outcome of circle Id looks difficult to determine which grid circles that the number of points (x,y,z) lies on it. I managed to extract the number of points (x,y,z) in one grid circle. The coding is shown below. However i have 7125 grid circles where I need to determine the number of points that lie on each grid circle. do you have any ideas?

% Input Number of Points
numPoints = 5475;
data = load('Sample.txt');
x = data (:,2);% x = lon
y = data (:,1);% y = lat
z = data (:,3);% z = H
plot(x, y, '.r');

% Input of Grids. Example of 1 point grid
x0 = 100.6;
y0 = 9.825;
R = 0.09; %Radius

% Plot circle
pos = [x0-R, y0-R, 2*R, 2*R];
rectangle('Position', pos, 'Curvature',[1 1]);
hold on;

% Determine how many points are in the circle .
count = sum(((x-x0).^2+(y-y0).^2<=R^2));
listC = count.';

TotalC = [y0 x0 listC];

%Extract Data in Circle
distances = sqrt((x-x0) .^ 2 + (y-y0) .^ 2);
InsideCircle = distances <=R;

xS = x(InsideCircle); % Selected Lon
yS = y(InsideCircle); % Selected Lat
zS = z(InsideCircle); % Selected H

%Save Total Count with Zero
fid = fopen('TotalCount.txt','wt');
fprintf(fid,'%.6f \t %.6f \t %.0f\n',[y0 x0 listC].');
fclose(fid);

% %Save Data Extraction Points in Circle Grid Point
append1 = [yS xS zS];
fid = fopen('1gridpts.txt','wt');
fprintf(fid,'%.6f \t %.6f \t %.5f\n',[yS xS zS].');
fclose(fid);

[hanif hamden]

i've modify some of your coding. How I want to declare the circleID belongs to which (x0,y0)?

[Cedric]

I don't fully understand what you are trying to achieve apparently. It seems to me that the last export that you are performing is quite similar to the one that I implemented, but that I implemented it already for managing all points and circles.

If you just want to save the counts in addition, here is an updated version that saves them in a second spreadsheet in the same workbook:

clc;clear;

% Load points.
numPoints = 5475 ;
data = load( 'Sample.txt' ) ;
x = data (:, 2 );
y = data (:, 1) ;
z = data (:, 3) ;                      
plot( x, y, '.r' ) ;

% Load circles.
grid = load( 'Grid025.asc' ) ;
[nCircles, ~] = size( grid ) ;
x0 = grid (:, 2) ;
y0 = grid (:, 1) ;
R = 0.09 ;

% Prealloc array of counts, and cell array for storing relevant point parameters per circle.
countPerCircle = zeros( nCircles, 1 ) ;
pointsPerCircle = cell( nCircles, 1 ) ;

% Initialize figure.
figure() ;
hold on ;

% Process circles.
for circleId = 1 : nCircles
    % Plot circle.
    pos = [x0(circleId)-R, y0(circleId)-R, 2*R, 2*R] ;
    rectangle( 'Position', pos, 'Curvature', [1, 1] ) ;
    
    % Compute distance from all points to current circle center.
    distances = sqrt( (x - x0(circleId)).^2 + (y - y0(circleId)).^2 ) ;
    
    % Flag and count relevant points.
    isInsideCircle = distances <= R ;
    countPerCircle(circleId) = sum( isInsideCircle ) ;
    
    % Extract relevant parameters.
    pointsId = find( isInsideCircle ) ;
    pointsX = x(isInsideCircle) ; % Selected Lon
    pointsY = y(isInsideCircle) ; % Selected Lat
    pointsZ = z(isInsideCircle) ; % Selected h
    
    % Store in relevant output cell. We build the first column as a vector of elements equal
    % to circleId, whose size matches the size of the other columns (e.g. pointsId).
    pointsPerCircle{circleId} = [circleId * ones( size( pointsId )), pointsId, pointsX, pointsY, pointsZ] ;
end

% Build data to export to Excel by concatenating all cells of the output cell array vertically
% into a large numeric array, converting it into a cell array, and concatenating to a first
% "header" row.
pointsPerCircle = num2cell( vertcat( pointsPerCircle{:} )) ;
pointsPerCircle = [{'CircleID', 'PointID', 'PointX', 'PointY', 'PointZ'}; pointsPerCircle] ;
 
% Export to Excel, spreadsheet "Points".
xlswrite( 'PerCircle.xlsx', pointsPerCircle, 'Points' ) ;
    
% Build output buffer for "count per circle".
countPerCircle = num2cell( [(1 : nCircles).', countPerCircle] ) ;
countPerCircle = [{'CircleID', 'Count'}; countPerCircle] ;

% Export to Excel, spreadsheet "Counts".
xlswrite( 'PerCircle.xlsx', countPerCircle, 'Counts' ) ;

This can be exported to text file if you prefer.

[hanif hamden]

It's okay, maybe I can modify some of the coding. One more question. How I would like to determine the circleID belongs to which grid points (x0,y0).

[Cedric]

Just update the code for building the output buffer for the counts, e.g. as follows:

% Build output buffer for "count per circle".
countPerCircle = num2cell( [(1 : nCircles).', x0, y0, countPerCircle] ) ;
countPerCircle = [{'CircleID', 'X0', 'Y0', 'Count'}; countPerCircle] ;

[hanif hamden]

for countPerCircle is okay because the circleID follows the nCircles. I can't solve for the output cell. I try to put like this..

output{circleId} = [circleId * ones( size( pointsId )), pointsId, x0(circleId), y0(circleId),  pointsX, pointsY, pointsZ]

This is the error..

Error using horzcat

Dimensions of matrices being concatenated are not consistent.

Error in DemoCircle2 (line 48)

output{circleId} = [circleId * ones( size( pointsId )), pointsId, x0(circleId), y0(circleId), pointsX, pointsY, pointsZ] ;

[Cedric]

It was a good attempt, but you tried to concatenate a scalar (e.g. x0(circleId) which is a single number) with column vectors, hence the error. If you look, circleId is also a scalar and I multiplied it by a vector of 1s whose size is compatible with other column vectors before concatenation. You have to do the same for x0 and y0.

I am updating the last version of the code, because it seems that you have been using the former one. Here, I create a column vector of 1s that I use for expanding circleId, X0, Y0:

clc ;
clear ;

% Load points.
numPoints = 5475 ;                  % *** Is this useful?
data = load( 'Sample.txt' ) ;
x = data (:, 2 ) ;
y = data (:, 1) ;
z = data (:, 3) ;                      
plot( x, y, '.r' ) ;

% Load circles.
grid = load( 'Grid025.asc' ) ;
[nCircles, ~] = size( grid ) ;
x0 = grid (:, 2) ;
y0 = grid (:, 1) ;
R = 0.09 ;

% Prealloc array of counts, and cell array for storing relevant point parameters per circle.
countPerCircle = zeros( nCircles, 1 ) ;
pointsPerCircle = cell( nCircles, 1 ) ;

% Initialize figure.
figure() ;
hold on ;

% Process circles.
for circleId = 1 : nCircles
    % Plot circle.
    pos = [x0(circleId)-R, y0(circleId)-R, 2*R, 2*R] ;
    rectangle( 'Position', pos, 'Curvature', [1, 1] ) ;
    
    % Compute distance from all points to current circle center.
    distances = sqrt( (x - x0(circleId)).^2 + (y - y0(circleId)).^2 ) ;
    
    % Flag and count relevant points.
    isInsideCircle = distances <= R ;
    countPerCircle(circleId) = sum( isInsideCircle ) ;
    
    % Extract relevant parameters.
    pointsId = find( isInsideCircle ) ;
    pointsX = x(isInsideCircle) ; % Selected Lon
    pointsY = y(isInsideCircle) ; % Selected Lat
    pointsZ = z(isInsideCircle) ; % Selected h
    
    % Build vector of 1s for expanding scalar to column vectors whose size
    % is compatible with other column vectors. Expand circleId, and relevant
    % elements of x0 and y0.
    onesVector = ones( size( pointsId )) ;
    circleIdExpanded = circleId * onesVector ;
    x0Expanded = x0(circleId) * onesVector ;
    y0Expanded = y0(circleId) * onesVector ;
    
    % Store in relevant output cell. We build the first column as a vector of elements equal
    % to circleId, whose size matches the size of the other columns (e.g. pointsId).
    pointsPerCircle{circleId} = [circleIdExpanded, x0Expanded, y0Expanded, pointsId, pointsX, pointsY, pointsZ] ;
end

% Build data to export to Excel by concatenating all cells of the output cell array vertically
% into a large numeric array, converting it into a cell array, and concatenating to a first
% "header" row.
pointsPerCircle = num2cell( vertcat( pointsPerCircle{:} )) ;
pointsPerCircle = [{'CircleID', 'X0', 'Y0', 'PointID', 'X', 'Y', 'Z'}; pointsPerCircle] ;
 
% Export to Excel, spreadsheet "Points".
xlswrite( 'PerCircle.xlsx', pointsPerCircle, 'Points' ) ;
    
% Build output buffer for "count per circle".
countPerCircle = num2cell( [(1 : nCircles).', x0, y0, countPerCircle] ) ;
countPerCircle = [{'CircleID', 'X0', 'Y0', 'Count'}; countPerCircle] ;

% Export to Excel, spreadsheet "Counts".
xlswrite( 'PerCircle.xlsx', countPerCircle, 'Counts' ) ;

[hanif hamden]

Now I get it.Thank you very much Sir. :)

[Cedric]

My pleasure! :-)