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

Question

Mini Me 2021-10-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1564221-number-of-points-inside-polygon-is-returning-0-but-should-be-a-bigger-number

评论： Mini Me 2021-10-27

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 个评论
显示 2更早的评论隐藏 2更早的评论

Mini Me 2021-10-15

编辑：Mini Me 2021-10-19

在 MATLAB Online 中打开

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 2021-10-19

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 2021-10-19

Files.zip

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 2021-10-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1564221-number-of-points-inside-polygon-is-returning-0-but-should-be-a-bigger-number#answer_811483

编辑：Cris LaPierre 2021-10-19

在 MATLAB Online 中打开

Files.zip

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

% 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 个评论
显示 4更早的评论隐藏 4更早的评论

Cris LaPierre 2021-10-27

在 MATLAB Online 中打开

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 2021-10-27

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.

请先登录，再进行评论。

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论