How to count pixel which are coming on the edge of the polygon?

Hello community,

I am facing issue in masking the data according to shape file. I took Indian latitude and longitude information. I have masked using the inpolygon function. Each pixel has size 2.5 * 2. I have a polygon, I want to select each pixel that placed in or on my polygon.The issue is that the pixel which is on the polygon is not considered by inpolygon command, however from the plot (left panel: after mask, right panel:after mask), one can readily see some pixel are cutting the polygon. I have shared the data and polygon info in the attachment.

I want each pixel which inside or atleast cutting the ploygon.

I would really appreciate any sort of help.

Thanking you in advance.

Answer 1

[KSSV]

Your data from .mat file already has lot of NaN values. Thee is no data in the mat file. Check your data.

The same code for some dummy data works fine.

lat1=0:2:40; 
lon=60:2.5:100;
% data = load ('data.mat') ;
shape=xlsread('indiaonlylatlon.xls');
[x,y]=meshgrid(lon,lat1);
[isin,ison] =inpolygon(y,x,shape(:,2),shape(:,1));
id=isin|ison;
% var1=data.ans;
var1=rand(size(x));
var1(~id)=NaN;
pcolor(x,y,var1)
hold on
plot(shape(:,1),shape(:,2),'r')

[Subhodh Sharma]

@KSSV Thanks for your response. I have updated the code, I hope this will make more sense to my need.

Yes the original data has lot of NaNs, because the observations are at scattered loaction. I can't do anything to it.

I want all the pixels which are cutting the edge of Indian boundary. The size of each pixel is 2.5 * 2. Maybe the 'inpolygon' considering only the center of pixel for calculating both 'isin' or 'ison' that's why if any pixels centre is going out the polygon, it's considering outside point.

I want to extract all those pixels whose one side or tiny portion is cutting the polygon.

I hope this makes more clear.

Thanking you once again for your consideration.

[KSSV]

That case go for knnsearch.

[KSSV]

lat=0:2:40; 
lon=60:2.5:100;

[nodes,X,Y,P] = Mesh(lon,lat) ;

var1 = load ('data.mat') ;
var1 = var1.ans ; 

shape=xlsread('indiaonlylatlon.xls');
[x,y]=meshgrid(lon,lat);
nel = size(nodes,1) ; 
iwant = zeros(nel,1) ;
for i = 1:size(nodes,1)
    [isin,ison] =inpolygon(shape(:,1),shape(:,2),X(:,i),Y(:,i));
    if any(isin|ison)
        iwant(i) = 1 ;
    end
end
idx = logical(iwant) ;

figure
hold on
plot(shape(:,1),shape(:,2),'r')
patch('faces',nodes(idx,:),'vertices',P,'facevertexcdata',var1(:),'facecolor','interp','edgecolor','k') ;


function [nodes,X,Y,P] = Mesh(lon,lat) 
Nx = length(lon)-1 ;   
Ny = length(lat)-1 ;       
% From here dont change
nel = Nx*Ny ;        % Total Number of Elements in the Mesh
nnel = 4 ;           % Number of nodes per Element
% Number of points on the Length and Breadth
npx = Nx+1 ;
npy = Ny+1 ;
nnode = npx*npy ;      % Total Number of Nodes in the Mesh
% Discretizing the Length and Breadth of the plate
[xx, yy] = meshgrid(lon,lat) ;
% To get the Nodal Connectivity Matrix
coordinates = [xx(:) yy(:)] ;
NodeNo = 1:nnode ;
nodes = zeros(nel,nnel) ;
% If elements along the X-axes and Y-axes are equal
if npx==npy
    NodeNo = reshape(NodeNo,npx,npy);
    nodes(:,1) = reshape(NodeNo(1:npx-1,1:npy-1),nel,1);
    nodes(:,2) = reshape(NodeNo(2:npx,1:npy-1),nel,1);
    nodes(:,3) = reshape(NodeNo(2:npx,2:npy),nel,1);
    nodes(:,4) = reshape(NodeNo(1:npx-1,2:npy),nel,1);
% If the elements along the axes are different
else%if npx>npy
    NodeNo = reshape(NodeNo,npy,npx);
    nodes(:,1) = reshape(NodeNo(1:npy-1,1:npx-1),nel,1);
    nodes(:,2) = reshape(NodeNo(2:npy,1:npx-1),nel,1);
    nodes(:,3) = reshape(NodeNo(2:npy,2:npx),nel,1);
    nodes(:,4) = reshape(NodeNo(1:npy-1,2:npx),nel,1);
end
P = [xx(:) yy(:)] ; 
%
% Plotting the Finite Element Mesh
% Initialization of the required matrices
X = zeros(nnel,nel) ;
Y = zeros(nnel,nel) ;
% Extract X,Y coordinates for the (iel)-th element
  for iel = 1:nel
      X(:,iel) = coordinates(nodes(iel,:),1) ;
      Y(:,iel) = coordinates(nodes(iel,:),2) ;
  end
end

And this will be the output of the code.