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

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

[KSSV]

We cennot check your code as the values are not given. Show us complete code with data, so that we can check.

[Mini Me]

The shapefile were uploaded as .txt file and the other .mat files Please change the .txt extension to .shp.

%% Loading the coordinates file
clear all

% 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'});
pfd_shp = shaperead("Block_9539_PFD_Contour.shp");

%% 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);
hold off

%% Plot Census bloks, pfd contour and population distribution
% 31 is block 9900, 42 is block 9539, 41 is block 9540
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');
hold on
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);

hold on

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(X_pop_trim,Y_pop_trim,Contour(l).latitude,Contour(l).longitude);
   
    pop_covered = numel(X_pop_trim(in));
   
   % disp(pop_covered)
   
    total_pop_covered = total_pop_covered + pop_covered;

end 

geobasemap('satellite')

hold off 

fprintf ('The total population covered in %s is %d people',Block_Name,total_pop_covered);

[Cris LaPierre]

I got the following errors:

Error using load

Unable to find file or directory 'Guam_PFD_Contour_1arc_noBerm262mtest.mat'.

Then

Error using openShapeFiles>checkSHP (line 91)

Invalid shapefile C:\Users\...\MATLAB Answers\tl_2019_66_tract.shp.

Error in openShapeFiles (line 19)

[basename, ext] = checkSHP(basename,shapeExtensionProvided);

Error in shaperead (line 210)

= openShapeFiles(filename,'shaperead');

Perhaps try zipping your 3 shapefiles and attaching that?

[Mini Me]

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

Answer 1

[Cris LaPierre]

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

[Mini Me]

Thanks for your update. My X_pop_trim had lon values in it and Y_pop_trim has lat values because inpolygon takes long before lats. Their name did not match the data inside. That was done on purpose. That's weird tahgt it works when you flip them.

What if I also wanted to get the area of each contour, would the polyarea function work?

[Cris LaPierre]

I agree. X is lon and Y is lat. inpolygon doesn't know what is lat and what is lon. What is important is that xq and yq be consistent with what your xv and yv inputs are. You use Contour(l).latitude,Contour(l).longitude for xv,yv so xq,yq needs to be (lat,lon), too.

• [in,on] = inpolygon(xq,yq,xv,yv)

I tried to create a geoscatter of your X_pop_trim and Y_pop_trim and got an error about lat values having to be less than 90. The X values are in the 140s. That is why I flipped them, and your code suddenly started to work.

I modified the code to color the circles inside the polygon green. It correctly identifies the points inside a polygon.

polyarea(x,y) returns the area of the 2-D polygon defined by the vertices in vectors x and y.

[Mini Me]

Ah I see your point, I never changed the contour(l).lat/lon to match X and Y_pop.

Thanks so much!

[Mini Me]

I used this line to get the area of the polygon contour but the answer does not seem accurate. 7.123e-8

Contour_area = polyarea (Contour(l).latitude,Contour(l).longitude); 

[Cris LaPierre]

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)

[Mini Me]

Interestingly enough I used this function "Geopolyarea" found on the file exchange to compare the results. https://www.mathworks.com/matlabcentral/fileexchange/75162-geopolyarea?s_tid=srchtitle_polyarea%20_3

I don't get the same results.