Select points inside a polygon
17 次查看(过去 30 天)
显示 更早的评论
I have 96 points (longitude and latitude); each point is a center of 0.5 x 0.5 pixel (box), on the other hand, I have a polygon, I want to select each pixel that placed in or on my polygon.
I can do it easily like this scrip below but the problem of this code is it just considers a point (center of the pixel) not an entire pixel to index.
polygon1_x = polygon1_x.'; % x of polygons
polygon1_y = polygon1_y.'; % y of polygons
lat = Points.lat; % x of points = lat
lon = Points.lon; % y of points = lon
[in,on] = inpolygon(lat,lon,polygon1_x,polygon1_y); % Logical Matrix
inon = in | on; % Combine ‘in’ And ‘on’
idx = find(inon(:)); % Linear Indices Of ‘inon’ Points
latcoord = lat(idx); % X-Coordinates Of ‘inon’ Points
loncoord = lon(idx); % Y-Coordinates Of ‘inon’ Points
clf
figure(1)
plot(lon, lat, '.') % Plot All Points
hold on
plot(polygon1_y, polygon1_x, '.') % Plot Polygon
plot(loncoord, latcoord, 'gp')
OUTPUT = idx; % the output is idx
So I would be grateful if anyone can told me how I can index the row number of points from Point.mat which placed in/on my shapefile if the points are center of 0.5 x0.5 pixels. In this way, I can effectively select pixels that are in or on my polygon. I need the output like idx.
Thank you so much.
采纳的回答
Bruno Luong
2020-8-4
编辑:Bruno Luong
2020-8-4
Here is the code using POLYSHAPE
load('Points.mat')
load('polygon1_x.mat')
load('polygon1_y.mat')
lat = Points.lat; % x of points = lat
lon = Points.lon; % y of points = lon
wpixel = 0.5;
hpixel = 0.5;
pixel0 = [-1 1 1 -1;
-1 -1 1 1]' .* [wpixel, hpixel]/2;
Island = polyshape([polygon1_y(:) polygon1_x(:)]);
for k=1:length(lat)
poly = polyshape([lon(k) lat(k)]+pixel0);
Ik = intersect(poly,Island);
isin = ~isempty(Ik.Vertices);
pixel(k).poly = poly;
pixel(k).isin = isin;
end
% Index of pixel that intersect with Island
idxin = find([pixel.isin]) % <=== HERE IS THE INDEX YOU NEED
% Graphical output
close all
plot(Island)
hold on
for k=1:length(pixel)
if pixel(k).isin
color = 'r';
else
color = 'w';
end
plot(pixel(k).poly, 'FaceColor', color);
end
axis equal
更多回答(2 个)
Bruno Luong
2020-8-4
But you already have the index in your code. here
...
inon = in | on; % Combine .in. And .on’
idx = find(inon(:));
...
What you think the idx are? Just use it
T(idx)
Precipitation(idx)
KSSV
2020-8-4
You can use the previous questions answer..with a slight change while saving the points.
clc; clear all ;
load("polygon1_x.mat") ;
load("polygon1_y.mat") ;
load("lon.mat") ;
load("lat.mat") ;
% remove Nan's from the data
xv = polygon1_y ; yv = polygon1_x ;
xv(isnan(xv)) = [] ;
yv(isnan(yv)) = [] ;
%
loncoord = cell([],1) ;
latcoord = cell([],1) ;
count = 0 ;
for i = 1:length(lon)
i
% make grid around (lon,lat)
x = lon(i)-0.5:0.5:lon(i)+0.5 ;
y = lat(i)-0.5:0.5:lat(i)+0.5 ;
[X,Y] = meshgrid(x,y) ;
[in,on] = inpolygon(X(:),Y(:),xv,yv); % Logical Matrix
inon = in | on; % Combine ‘in’ And ‘on’
idx = find(inon(:)); % Linear Indices Of ‘inon’ Points
if any(idx)
count = count+1 ;
loncoord{count} = X(idx); % Y-Coordinates Of ‘inon’ Point
latcoord{count} = Y(idx) ;
end
end
plot(polygon1_y,polygon1_x,'b')
hold on
for i = 1:length(loncoord)
plot(loncoord{i},latcoord{i},'*r')
end
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Polygons 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!