4

I use Mapbox to display some moving markers over a map with several polygons drew over it (feature layers).

I have serveral markers (~1000) and multiple feature layers (polygons ~200) that represent different regions. The markers are moving (once per second with a constant velocity) and I have to know at any moment on which polygon does the marker position reside on.

First idea is the one to iterate for each marker every second over all polygons and check with a point in polygon formula (for geopositions?) if it is inside or not. This solution will not work as the complexity of O(NrMarkers * NrPolygons) seems too mcuh to run every second (might work but will put too mcuh stress on the server that moves the marker positions).

Now, I think maybe there is some way I could use more complex structures like k-d trees but it will be pretty hard to implement them and still there would be pretty much computation to be done each second.

I was also thinking of pre-processing because the polygons are known from the start and never change. Maybe I could somehow store for some known coordinates what polygons intersect those coordinates and use this info to exclude some polygons (or directly know the correct one) for any given position query.

So, to formulate it as an algorithmic question:

Given a list of P polygons (~200) (each having between 3 and 20 vertexes) and a list of N points (~1000) (which can either be still or have a constant speed that might change direction from time to time) find an algorithm that tells for each of the N points with which polygon shape does it intersect (inside which polygon is it, if any). This query will be asked every second for the same list of polygons and the same list of points (points which have their position updated after each second).

11
  • I'd be looking at GIS extensions to existing databases. For example, postgresql and this query.
    – user40980
    Commented Oct 13, 2015 at 18:11
  • So are the polygons guaranteed to no overlap?
    – MetaFight
    Commented Oct 13, 2015 at 18:15
  • @MetaFight They are not at the moment but if this simplifies the solution I can use only polygons that don't overlap, I can make my app work without overlapping areas.
    – Cristy
    Commented Oct 13, 2015 at 19:51
  • 2
    @Cristy a database query with geospatial indexes can be very fast. And its a solved problem meaning you don't have to spend any more time trying to do it yourself.
    – user40980
    Commented Oct 13, 2015 at 19:59
  • 1
    I'm not a noder, so I'm not completely familiar with how to hook up node to postgresql. I do hear that noders have a particular affinity for Mongo, so you might try Geospatial Indexes and Queries from Mongo to see if that does what you want. You could possibly do something with couch too. Glance at GIS.SE's nosql tag also.
    – user40980
    Commented Oct 13, 2015 at 20:04

2 Answers 2

2

This geographical problem is called Point in Polygon, and there's a solution in Python here.

For the sake of completeness here it is in JavaScript with some examples:

// x,y are point coordinates
// poly is a list of tuples [[x,y], [x,y], ....]
function isPointInPolygon(x,y,poly) {
    var num = poly.length,
        i = 0;
        j = num - 1;
        c = false;

    for( i = 0; i < num; i++ ) {
        if( ((poly[i][1] > y) != (poly[j][1] > y)) && (x < (poly[j][0] - poly[i][0]) * (y - poly[i][1]) / (poly[j][1] - poly[i][1]) + poly[i][0]) ) {
            c = !c;
        }
        j = i;
    }
    return(c);
}

isPointInPolygon(5,1, [[0,0], [7,3], [9,0]]); // true
isPointInPolygon(5,3, [[0,0], [7,3], [9,0]]); // false

This shouldn't be used for lat/lng coordinates for very large areas. City scale should be fine though.

Update

Rereading the question I'm not sure I've sufficiently answered. I've added new code where the client's browser will run this much better than O(NrMarkers * NrPolygons) as long as your polygons don't overlap too much. This expanded code creates, and can store, bounding rectangles for polygons to reduce expensive calls to isPointInPolygon.

// Not to be used for large areas or near equator or 0/180 longitude
// x,y are point coordinates
// poly is a list of tuples [[x,y], [x,y], ....]
function isPointInPolygon(x, y, poly) {
  var num = poly.length,
      i = 0;
  j = num - 1;
  c = false;

  for (i = 0; i < num; i++) {
    if (((poly[i][1] > y) != (poly[j][1] > y)) && (x < (poly[j][0] - poly[i][0]) * (y - poly[i][1]) / (poly[j][1] - poly[i][1]) + poly[i][0])) {
      c = !c;
    }
    j = i;
  }
  return (c);
}

function arePointsInPolygons(points, polygons) {
    var pointnum = points.length,
        polynum = polygons.length,
        vertnum,
        pb = polybounds,
        results = [],
        pt, pg, v;

    if( pb.length === 0 ) {
        // get the polygon bounds for each polygon
        for( pg = 0; pg < polynum; pg++ ) {
            pb.push([false, false, false, false]);
            vertnum = polygons[pg].length;

            for( v = 0; v < vertnum; v++ ) {
                // north
                if( pb[pg][0] === false || polygons[pg][v][1] > pb[pg][0] ) {
                    pb[pg][0] = polygons[pg][v][1];
                }
                // east
                if( pb[pg][1] === false || polygons[pg][v][0] > pb[pg][1] ) {
                    pb[pg][1] = polygons[pg][v][0];
                }
                // south
                if( pb[pg][2] === false || polygons[pg][v][1] < pb[pg][2] ) {
                    pb[pg][2] = polygons[pg][v][1];
                }
                // west
                if( pb[pg][3] === false || polygons[pg][v][0] < pb[pg][3] ) {
                    pb[pg][3] = polygons[pg][v][0];
                }
            }
        }
        polybounds = pb;
    }

    // check which points fall within the bounds, if yes call isPointInPolygon
    for( pt = 0; pt < pointnum; pt++ ) {
        results.push([]);
        for( pg = 0; pg < polynum; pg++ ) {
            if( points[pt][0] > pb[pg][3] && points[pt][0] < pb[pg][1] &&
               points[pt][1] > pb[pg][2] && points[pt][1] < pb[pg][0] ) {
                if( isPointInPolygon(points[pt][0], points[pt][1], polygons[pg]) ) {
                    results[pt].push(pg);
                }
            }
        }
    }

    return(results);
}

// Test
var i,
    inside,
    numiter = 1000000,
    points = [];
for( i = 0; i < numiter; i++ ) {
    points.push([parseInt(Math.random()*10,10), parseInt(Math.random()*10,10)]);
}

// test isPointInPolygon for list of many points
inside = 0;
var start = new Date().getTime();
for( i = 0; i < numiter; i++ ) {
    inside += isPointInPolygon(points[i][0], points[i][1], [[0, 0], [7, 3], [9, 0]]) ? 1 : 0;
}
console.log(new Date().getTime() - start);


// do same test with optimization
var polybounds = [];
var start = new Date().getTime();
arePointsInPolygons(points, [[[0, 1], [7, 3], [9, 0]],[[0, 3], [3, 3], [3, -1]]]);
console.log(new Date().getTime() - start); // slower, builds polybounds

var start = new Date().getTime();
arePointsInPolygons(points, [[[0, 1], [7, 3], [9, 0]],[[0, 3], [3, 3], [3, -1]]]);
console.log(new Date().getTime() - start);

// arePointsInPolygons returns array of length equal to number of points containing arrays of polygon indexes in which the points fall - empty indicates not in any.

Limitations update

Looking at the map below: USA at different projections (source)

The northern border, which follows the 49th parallel, is sometimes straight and sometimes curved on the maps. For web maps the projection used shows the border as horizontal. The danger lies when a point falls within the area in one projection but not in the other. The path along a curved surface is ambiguous. See diagram:

Polygons contain different space depending on projection

5
  • How does this answer my question? I didn't ask how to check if a point resides in a polygon, the question is about a good complexity algorithm that does what is explained in the question. "This shouldn't be used for lat/lng coordinates for very large areas. City scale should be fine though." How does the range of the coordinates increase the complexity, I don't get it.... ?
    – Cristy
    Commented Oct 14, 2015 at 9:41
  • Yeah sorry, I wrote it late last night. I was just updating the answer. Because the earth is round, lines between two points is ambiguous. Imagining a situation where you provide a large rectangular polygon, this provided method will determine if points fall within the top and bottom latitudes. So on the globe it will check a shape similar to a windshield rather than a trapezoid which may be expected. I'll update the answer with an image.
    – Cyrille
    Commented Oct 14, 2015 at 10:55
  • Thanks for the answer update, but as I said, I don't want a O(NrMarkers * NrPolygons) solution as it is too slow to be done in real-time each second. As for the check precision, it doesn't have to be very exact (eg: it doesn't matter that the earth is actually a sphere)
    – Cristy
    Commented Oct 14, 2015 at 16:21
  • @Cristy unless you specify a lot more information about the nature of the points and polygons that will let us get a better idea of the relationships between each O(points*polys) is the best you will do. As you say that the location of the polygons can change from one time to the next, it is likely that O(points*polys) is the best that you will ever be able to do without going through and implementing your own GIS enabled database indexes - not an easy task.
    – user40980
    Commented Oct 14, 2015 at 21:01
  • @MichaelT Well, I didn't ask for a one line solution, I said that even implementing my own k-d trees might be a solution if it is optimal or improves the running time a lot. And I can't accept this as the best solution as it doesn't really answer anything that isn't already mentioned or described in the question. I would rather accept your comment mentioning the geospatial query in MongoDB.
    – Cristy
    Commented Oct 15, 2015 at 15:17
1

Hopefully not too late... I'm having a similar problem, except that the number of my points are in billions, and the number of polygons could be in thousands.

While this idea did not help me, it might help you: if the area covered by the polygons is relatively small and/or the resolution is not too fine, you can create a grid of all points in the area, and for each point pre-compute the polygon it belongs in. Then, at run time, you'd find the 4 grid points that span a square/rectangle that contains your query point, and for each grid point you'll figure out its polygon, then you'll have 4 polygons to work with.

On the whole, I agree with @MichaelIT that spatial databases are pretty fast, and PostGIS (PostgreSQL extension) can easily cache in memory the shapes of polygons you work with, so, it won't be any slower than your custom code.

1
  • I've ended up using the GeoSpatial database queries from MongoDB and for the number of polygons I have it works fast enough (completes in the same time as a simple select query does).
    – Cristy
    Commented Dec 3, 2015 at 13:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.