How to calculate the points of a line which lies inside polygons?

Hi all,

I have GPS coordinates of an off-road vehicle and I have to check when the vehicle lies on a road or on a field. I have downloaded the shape file of my region, where there are the bounding box of each field. Below, you will find my code, that is very slow because in the shape file there are around 3e5 polygons. How could I speed up the code? In my opinion, the calculation could be speeded up, if I can avoud the inner for, but I do not know how can I do it.

for i=1:length(GPSData)
   OnField{i}=zeros(size(GPSData(i).Longitude));         
    for iS=1:length(FieldUseSHP)
        a=inpolygon(GPSData(i).Longitude*pi/180,GPSData(i).Latitude*pi/180,FieldUseSHP(iS).Lon,FieldUseSHP(iS).Lat);
        if any(a)==1
            OnField{i}=iS.*a;
        end
    end
end

Thank you in advance

Best regards,

Pietro

[Adam Danz]

If the road has a fixed width, and if you have the coordinates for the center of the road, this problem becomes much easier.

In that case, you can just calculate the lateral deviation between the vehicle and the road's center and if it exceeds half of the road width, you know the vehicle is off of the road.

So 1) if the road width is fixed and known and 2) you have the center coordinates of the road, I'd be glad to walk you through the rest of that.

[pietro]

Hi Adam,

yes, I have the coordinates of the road and of the field profiles. The problem is solved as you suggested, but the problem remains with the field. Here a picture of how my problem looks like:

https://i.postimg.cc/Kz7JqMCM/untitled.jpg

as you can see the tractor line is enclosed by a field profile. In that case, I have used the approach shown above, which is very slow.

Pietro

[Adam Danz]

I see. What does the data look like that produces the image you shared? In other words, the brown field contains ~20 line segments and the yellow field contains 4 segments. Do you have the 20 (x,y) coordinates that produce the brown field and the 4 (x,y) coordinates that produce the yellow field?

If yes, are those coordinates used to define the polygon used in inpolygon()?

Also, what are your two loops? You should have a set of (x,y) coordinates that define the tractor's trajectory and those whould be the first inuts to the inpolygon() function.

[pietro]

Hi,

I have several GPS vehicle data, each is an acquisition and the first for loops through all the acquisitions. The other for loops through all the polygons of field plot contours.

In the plot I have shown only two field countours, that are the closest to the chosen acquisition.

Cheers

Pietro

Answer 1

[Matt J]

fconv=@(z) z(:).';
Longitude=cell2mat( cellfun( fconv, {GPSData.Longitude},'uni',0)  )*pi/180; 
Latitude=cell2mat( cellfun( fconv, {GPSData.Latitude},'uni',0)  )*pi/180; 

Onfield=cell(1,length(FieldUseSHP));

    for iS=1:length(FieldUseSHP)
        
        a=inpolygon(Longitude, Latitude,FieldUseSHP(iS).Lon,FieldUseSHP(iS).Lat);
 
        OnField{iS}=iS.*a;

    end