Number of points inside polygon is returning 0 but should be a bigger number

3 views (last 30 days)
Hello, I am working with this code where I generate random points inside a polygon then I superpose another polygon inside the first one and check for number of points that fall inside a the second polygon. The number of points return 0 when it should not be. The data plots fine. but the result value is 0 for both in and on. The code below shows the loop that does that. X_pop_trim,Y_pop_trim already defined. See code below
  4 Comments
Cris LaPierre
Cris LaPierre on 19 Oct 2021
I got the following errors:
Error using load
Unable to find file or directory 'Guam_PFD_Contour_1arc_noBerm262mtest.mat'.
Then
Invalid shapefile C:\Users\...\MATLAB Answers\tl_2019_66_tract.shp.
[basename, ext] = checkSHP(basename,shapeExtensionProvided);
Error in shaperead (line 210)
= openShapeFiles(filename,'shaperead');
Perhaps try zipping your 3 shapefiles and attaching that?
Mini Me
Mini Me on 19 Oct 2021
Here are the files. load Guam_PFD_Contour_1arc_noBerm262mtest.mat will not work since that file is too big.
I made changes to the code to only to load up the Lat, lon and DeltaPi separately.
See code above

Sign in to comment.

Accepted Answer

Cris LaPierre
Cris LaPierre on 19 Oct 2021
Edited: Cris LaPierre on 19 Oct 2021
In line 114 of the code you have shared, you have reversed your lat (you name it Y) and lon (you name it X) data. Flip those, and it finds people.
I've cleaned up some of your hold on/offs and used the figure command to separate your cartesian and geoplots into separate figures.
% Just to make the files available here
unzip('Files.zip')
% Your code #####################
%% Loading the coordinates file
% City Name
City_Name = 'Guam';
Block_Name ='Block 9539';
% Easth station Location
ES_lat = 13.416944;
ES_lon = 144.751944;
% Load the file that contains the pfd contour coordinates
%load Guam_PFD_Contour_1arc_noBerm262mtest.mat;
load Lat_T.mat;
load Long_T.mat;
load DeltaPi.mat;
% Load the file that contains the Guam Census Block Contour
Guam_shp = shaperead('tl_2019_66_tract.shp','Attributes',{'NAME','LAT','LON'});
Warning: Failed to open file /users/mss.system.Znncyb/tl_2019_66_tract.shx or file /users/mss.system.Znncyb/tl_2019_66_tract.SHX. Will build index from SHP file.
Warning: Failed to open file /users/mss.system.Znncyb/tl_2019_66_tract.dbf or file /users/mss.system.Znncyb/tl_2019_66_tract.DBF. Shape output structure will have no attribute fields.
pfd_shp = shaperead("Block 9539 PFD Contour.shp");
Warning: Failed to open file /users/mss.system.Znncyb/Block 9539 PFD Contour.shx or file /users/mss.system.Znncyb/Block 9539 PFD Contour.SHX. Will build index from SHP file.
Warning: Failed to open file /users/mss.system.Znncyb/Block 9539 PFD Contour.dbf or file /users/mss.system.Znncyb/Block 9539 PFD Contour.DBF. Shape output structure will have no attribute fields.
%% Determine the contour for threshold 0
M=contour(long_T,lat_T,DeltaPi,[0,0]);
M=M';
%% Put that in a better format
nulls=find(M(:,1)==0);
for k=1:numel(nulls)-1
Contour(k).longitude=M(nulls(k)+1:nulls(k+1)-1,1);
Contour(k).latitude=M(nulls(k)+1:nulls(k+1)-1,2);
num_elem(k)=numel(Contour(k).longitude);
end
%% Eliminate the smaller contours (here lower than 10 points)
threshold=10;
Contour=Contour(num_elem>=threshold);
num_elem=num_elem(num_elem>=threshold);
%% Sort by number of points of the contour
[~,ind_num]=sort(num_elem,'descend');
Contour=Contour(ind_num);
%% Plot Census bloks, pfd contour and population distribution
% 31 is block 9900, 42 is block 9539, 41 is block 9540
figure
S = geoplot(Guam_shp(41).Y,Guam_shp(41).X,'DisplayName',Block_Name,"LineWidth",3,"LineStyle","-","Color",'r','Marker','none');
hold on
T = geoplot(Guam_shp(42).Y,Guam_shp(42).X,'DisplayName',Block_Name,"LineWidth",3,"LineStyle","-","Color",'r','Marker','none');
geoplot (ES_lat,ES_lon,'DisplayName','Guam ES Gateway','Color','r','Marker','*','MarkerSize',20)
%% Creating the random population grid inside the censis block
B_9540 = Guam_shp(41).BoundingBox;
B_9540_minlat = min(B_9540(:,2)); %One value
B_9540_maxlat = max(B_9540(:,2));% One Value
B_9540_latspan = B_9540_maxlat - B_9540_minlat;
B_9540_minlon = min(B_9540(:,1));
B_9540_maxlon = max(B_9540(:,1));
B_9540_lonspan = B_9540_maxlon - B_9540_minlon;
B_9540_X = Guam_shp(41).X; %range of Lat
B_9540_Y = Guam_shp(41).Y; %range of Lon
B_9540_population = 2257;
X_pop = zeros(1,B_9540_population);
Y_pop = zeros(1,B_9540_population);
for i=1:B_9540_population
flagIsIn = 0;
while ~flagIsIn
x(i) = B_9540_minlon + rand(1) * B_9540_lonspan;
y(i)= B_9540_minlat + rand(1) * B_9540_latspan;
flagIsIn = inpolygon(x(i),y(i),B_9540_X,B_9540_Y);
if flagIsIn ==1
X_pop = [X_pop x(i)];
Y_pop = [Y_pop y(i)];
end
end
end
% Removing zeros from Arrays
X_pop_trim = X_pop(X_pop~=0);
Y_pop_trim = Y_pop(Y_pop~=0);
geoscatter(Y_pop_trim,X_pop_trim,'color','r','Marker','o')
%% Add the pfd contour on population coverage
total_pop_covered = 0;
for l = 1: size (Contour,2)
geoplot(Contour(l).latitude,Contour(l).longitude,"LineWidth",1,"LineStyle","-","Color",'c','Marker','none')
[in,on] = inpolygon(Y_pop_trim,X_pop_trim,Contour(l).latitude,Contour(l).longitude);
%>>>> ^ ^ ################## Reverse X and Y
pop_covered = numel(X_pop_trim(in));
total_pop_covered = total_pop_covered + pop_covered;
end
hold off
geobasemap('satellite')
fprintf ('The total population covered in %s is %d people',Block_Name,total_pop_covered);
The total population covered in Block 9539 is 159 people
  6 Comments
Cris LaPierre
Cris LaPierre on 27 Oct 2021
I think there is a units issue. The better function is probably areaint.
  • area = areaint(lat,lon,units)
You can specify units using the earthRadius function.
Note: both of these functions are part fo the Mapping Toolbox
% To obtain ara in km^2
rKm = earthRadius('km');
Contour_area = areaint(Contour(l).latitude,Contour(l).longitude,rKm)

Sign in to comment.

More Answers (0)

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!