function IN = inside_shape_proj(E,N,shapefile)
%% FIND ALL POINTS WITHIN A POLYGON shapefile
% USAGE NOTES:
% - Shapefile must be in map projection (avoids dateline problems)
% - Takes long for large data sets, e.g. 100000 points ~3 min
% - Many small polygons run much faster than few large ones
% - Stepwise search size is 100x100 km by default
%% USER INPUT
% E = east or x coordinate in same unit and projection as shapefile
% N = north or y coordinate in same units and projection as shapefile
% shapefile = full path to polygon shapefiles
% (input as cells with each cell containing the shapefile path+filename)
% (NB! Strings need to be converted to cells, e.g. cellstr('filename'))
% example usage:
% E = [x1 x2];
% N = [y1 y2];
% shapefile = ...
% {'/Documents/Antarctica/moa/outlines/moa_iceshelves_pol.shp'};
% output:
% IN = 1 (inside polygons), IN = 0 (outside polygons)
% step size for distributed search
step = 100000; %unit meters, default 100 km
% find lat and lon limits of input points
E_lim = [min(E)-step max(E)+step];
N_lim = [min(N)-step max(N)+step];
% E and N need to be column vectors
if sum(size(E)>1)>1 || sum(size(N)>1)>1
error ('E and N need to be vectors')
elseif sum(sum(size(E) - size(N))) > 0
error ('E and N vectors must be the same size')
elseif size(N,2) > 1
E = E';
N = N';
FLIP = true;
else
FLIP = false;
end
% initialize
IN = zeros(size(E)) == 1;
tic
for sf = 1 : length(shapefile)
S = shaperead(shapefile{sf});
% check if single part shapefile
if length(S) > 1
% Convert polygon parts from cell arrays to vector form
[Y X] = polyjoin({S(:).Y},{S(:).X});
else
% extract E and N data
X = S.X;
Y = S.Y;
end
clear S
disp(['Determining points in shapefile ',num2str(sf),' of ',num2str(length(shapefile))])
N_step = floor(N_lim(1)) : step : ceil(N_lim(2));
IN1 = zeros(length(N),length(N_step)) == 1;
% loop for each step length in north direction
for N_id = 1 : length(N_step)
N0 = N_step(N_id);
% this script can take a while, uncomment the line below to tack
disp([' Searching data in segment ',num2str(N_id),' of ',num2str(length(N_step)),''])
Epol = [E_lim(1) E_lim(2) E_lim(2) E_lim(1) E_lim(1)]; %E-coord for search rectangle
Npol = [N0+step N0+step N0 N0 N0+step]; %N-coord for search rectangle
inNcol = inpolygon(E,N,Epol,Npol);
% reduce data in loop and clip to north bounds
X0 = E(inNcol);
Y0 = N(inNcol);
IN0 = zeros(sum(inNcol),1) == 1;
% extract common polygon
[E_in1,N_in1] = polybool('intersection',Epol,Npol,X,Y);
% loop for each step size in east direction
for E0 = floor(E_lim(1)) : step : ceil(E_lim(2))
Epol = [E0 E0+step E0+step E0 E0];
inCell = sum(inpolygon(X0,Y0,Epol,Npol));
if inCell == 0
% if no points in rectangle then jump to next iteration
continue
else
% extract common polygon
[E_in,N_in] = polybool('intersection',Epol,Npol,E_in1,N_in1);
if isempty(E_in)
continue
else
IN0 = IN0 | inpolygon(X0,Y0,E_in,N_in);
end
end
end
IN1(inNcol,N_id) = IN0;
end
IN = sum(IN1,2) | IN;
end
if FLIP
IN = IN';
end