How do I get specific data using inpolygon in Matlab?
7 views (last 30 days)
Show older comments
harman bhaveja
on 8 Apr 2019
Commented: harman bhaveja
on 20 Apr 2019
So,
I've written a code where I'm reading flight data from an excel sheet and trying to find all the flights within a particular bound using the inpolygon code. I to only record want the flight data which has all its cordinates within the polygons bound, and if it is not within the bound I don't want it to be recorded, but it isn't working...
This is the code:
clc;
clear all;
if exist('coords.mat','file')
load coords
else
lat = []; % latitude
lon = []; % longitude
alt = []; % altitude
spd = []; % speed
flg = {}; % flight
tim = []; % time
end
lat1 = xlsread('plot.xlsx', 'C2:C86400');
lon1 = xlsread('plot.xlsx', 'D2:D86400');
alt1 = xlsread('plot.xlsx', 'E2:E86400');
flg1 = xlsread('plot.xlsx', 'A2:A86400');
lat10 = [6.6038,21.841667]; lon10 = [94.4166,62.375000];%P574
% Save and continue
save('coords.mat','lat','lon','alt','flg')
pause(1)
%% Render the recorded data
% Prepare figure
d=37040;
load coords %#ok<*UNRCH>
close all
% Settings
centerLoc = [12.9716,77.5946]; % LEVC
% Prepare UTM scenario
mstruct = defaultm('utm');
mstruct.zone = utmzone(centerLoc(1),centerLoc(2));
mstruct = defaultm(mstruct);
% Plot land contours
SHPdir = '.\SHPs\';
countries = shaperead([SHPdir 'ne_10m_admin_0_countries.shp'],...
'Selector',{@(x) strcmpi(x,'es'),'foo'},'UseGeoCoords', true);
% Change 'ES.VC' for the provinces/states of your preference or use a RegExp
% for all provinces: @(x) strcmpi(x,'ES.VC') => @(x) ~isempty(regexpi(x,'^ES.*$'))
provinces = shaperead([SHPdir 'ne_10m_admin_1_states_provinces.shp'],...
'Selector',{@(x) strcmpi(x,'ES.VC'),'region_cod'},'UseGeoCoords', true);
[x,y] = mfwdtran(mstruct,[countries.Lat provinces.Lat],[countries.Lon provinces.Lon]);
[xc,yc] = mfwdtran(mstruct,centerLoc(1),centerLoc(2));
[x1,y1]= mfwdtran(mstruct,lat1,lon1);
[x10, y10]= mfwdtran(mstruct,lat10,lon10);
X = [x10,fliplr(x10)-d];
Y = [y10,fliplr(y10)];
[in,on] = inpolygon(x1,y1,X,Y);
% (in) contains the indices of the points inside the parallelogram
% (on) contains the indices of the points on the parallelogram
j =1;
i =1;
p = 1;
flightno = [] ;
k = 0;
arr = 0;
for(flg1 = 1:86413)
while(i>1)
if(flg1(i)~=flg1(i-1))
j =j+1;
p = 1;
if (k ==1)
arr=arr+1;
%plot flight no. (j - 1)
flightno(arr) = (j-1);
k = 0;
end
end
end
%checking
%define polygon and write check condition
if(x1(i)~= x1(~in))
if (p == 1)
k = 1;
else
k = 0;
p = 0;
end
end
end
Can you help me out with this please?
0 Comments
Accepted Answer
Guillaume
on 8 Apr 2019
I don't have the mapping toolbox and knows nothing about it, so possibly the problem is with your coordinate transformation code.
With regards to finding which flight is entirely within a polygon here is how to go about it:
journeys = readtable('plot.xlsx'); %much simpler than multiple xlsread
%your mapping toolbox code here. I've no idea if it's correct
%...
[journeys.x, journeys.y] = mfwdtran(mstruct, journeys.latitude, journeys.longitude); %transform the flights latitude and longitude
in = inpolygon(journeys.x, journeys.y, X, Y);
[groupid, flights] = findgroups(journeys.flight); %assign unique id to each flight and apply to rows of the table
isallin = splitapply(@all, in, groupid); %are ALL points of the flight in the polygon? logical output
selectedflights = flights(isallin); %list of flights where all points are in the polygon
selectedjourneys = journeys(ismember(journey.flight, selectedflights), :); %portion of the table with only the selected flights.
9 Comments
Guillaume
on 18 Apr 2019
The correct line should have been:
inratio = splitapply(@(in) nnz(in)/numel(in), in, groupid); %ratio of points in polygon to points in flight
It's not my code that creates this empty table since my code does not create a m variable. All my variables have meaningful names so it's clear what their purpose is.
Anyway, the easiest way for you to find out what has gone wrong with your code is to debug it Step through the code one line at a time, see how the variables change and if they change the way you expected. If not, fix the line.
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!