46 views (last 30 days)

Show older comments

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
on 4 Aug 2020

Edited: Bruno Luong
on 4 Aug 2020

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

Bruno Luong
on 4 Aug 2020

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
on 4 Aug 2020

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

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

Start Hunting!