Here is the code using POLYSHAPE
load('Points.mat')
load('polygon1_x.mat')
load('polygon1_y.mat')
lat = Points.lat;
lon = 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
idxin = find([pixel.isin])
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
Best Answer