Number of points inside polygon is returning 0 but should be a bigger number
5 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
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
Invalid shapefile C:\Users\...\MATLAB Answers\tl_2019_66_tract.shp.
[basename, ext] = checkSHP(basename,shapeExtensionProvided);
= openShapeFiles(filename,'shaperead');
Perhaps try zipping your 3 shapefiles and attaching that?
採用された回答
Cris LaPierre
2021 年 10 月 19 日
編集済み: Cris LaPierre
2021 年 10 月 19 日
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'});
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);
%% 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);
6 件のコメント
Cris LaPierre
2021 年 10 月 27 日
- area = areaint(lat,lon,units)
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)
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Map Display についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!