[GIS] Checking if point is inside irregular polygon

javascriptpoint-in-polygonpolygon

I'm using Javascript (https://github.com/substack/point-in-polygon) to test if a point is inside polygon, however I think it's adapted from https://www.ecse.rpi.edu/~wrf/Research/Short_Notes/pnpoly.html

Following are the list of points,

let points = [
            [ 42.34624, -71.06024 ],
            [ 42.33558, -71.06616 ],
            [ 42.32632, -71.05835 ],
            [ 42.32987, -71.05428 ],
            [ 42.34732, -71.05432 ],
            [ 42.34618, -71.05973 ],
            [ 42.34624, -71.06024 ]
        ];
let testPoint = [
            [42.33288, -71.05804]

        ];`

let isInside = pp(testPoint, points);
expect(isInside).to.equal(true); // always fails

I'm not sure why it fails though. Is it because the algorithm doesn't cope with this "type" of polygon? Any hints please? Please point me to any alternatives to this approach as well.

Best Answer

If you are using Openlayers see here.
If you are using leaflet see here.
For pure javascript see here or for NPM resources here or here

There are other solutions beyond these but this should be enough to be getting on with.

Related Solutions

Polygon – Best Way to Determine if a Point is Within a Certain Distance of a GEOJSON Polygon

This really depends on a number of factors.

Are you talking about a large number of features? Do you have 100 neighborhoods and 1000 point locations, or a thousand times that number?
Is your precision requirement high enough that you'll need to take Great Circle calculations into account? At less than 10 miles the error will be around 10-30 meters, depending on your latitude (see 2).

If either of the above are true, then you'd be much better off parsing the GeoJSON into a spatial database and making the calculations from there. MongoDB offers some geospatial capability http://docs.mongodb.org/manual/core/2dsphere/ without having to mess around with GeoJSON parsing, but you'd have to play around with some of the spatial query operators to find something that suits. The other more established Spatial databases like PostGIS and Oracle will have more tools for the job (e.g. Oracle's SDO_WITHIN_DISTANCE), but you'll have more work getting up and running.

If you don't have loads of features and you're not worried about the great circle errors, then you could try using JSTS Topology Suite to help.

Here's a basic example. You'll need to work to calculate the basic lat/lng distances involved.

var tolerance = 10;

var geojsonPolygon = {
        "type": "Polygon",
        "coordinates": [
            [
                [10.0, 10],
                [20.0, 10],
                [20.0, 20],
                [10.0, 20],
                [10.0, 10]
            ]
        ]
};

var geojsonPoint =  {
        "type": "Point",
        "coordinates":
                [0, 0]
};

var reader = new jsts.io.GeoJSONReader();

var polygon = reader.read(geojsonPolygon);

var point = reader.read(geojsonPoint);

var ptDist = new jsts.algorithm.distance.PointPairDistance();
jsts.algorithm.distance.DistanceToPoint.computeDistance(polygon,point.getCoordinate(),ptDist);

if(ptDist.getDistance() < tolerance){
    // It's within 10
} else {
    // It's further than 10
};

Hope that helps.

[GIS] Point in Polygon – geojson – using Shapely / Python returning incorrect results

Shapely manual says "All operations are performed in the x-y plane."

in def getPOIs change

POIs.append((float(latitude), float(longitude)))

POIs.append((float(longitude), float(latitude)))

Here's the code that works; commented out all but one csv

import urllib2, json, csv
import numpy

from shapely.geometry import shape, Point


def readJson(url):
    """
    Returns a json file specified in @url.
    """
    response = urllib2.urlopen(url)
    return json.loads(response.read())

def readCSV(url):
    """
    Returns a csv file specified in @url.
    """
    response = urllib2.urlopen(url)
    return csv.DictReader(response, delimiter=',')


def getRegions():
    """
    Returns a dictionary formed by the id of a region and its coordinates.
    """
    dict = {}

    url = "https://nycdatastables.s3.amazonaws.com/2013-08-19T18:22:23.125Z/community-districts-polygon.geojson"
    data = readJson(url)
    for district in data['features']:
        dict[district['id']] = district['geometry']

    return dict


def getPOIs():
    """
    Returns a list of tuples of POIs lat/long coordinates.
    """
    urls = [
        "https://nycdatastables.s3.amazonaws.com/2013-06-04T18:02:56.019Z/museums-and-galleries-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-12-16T21:49:55.716Z/nyc-parking-facilities-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-20T16:06:05.136Z/mapped-in-ny-companies-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-11T18:59:27.269Z/nyc-public-school-locations-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-07-29T15:49:03.498Z/nyc-private-school-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-07-01T16:25:00.297Z/nyc-special-education-school-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-05T14:35:56.387Z/basic-description-of-colleges-and-universities-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-05T20:25:17.301Z/operating-sidewalk-cafes-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-04T14:40:48.764Z/community-health-centers-results.csv",
        #     "http://data.nycprepared.org/ar/dataset/dycd-after-school-programs-housing/resource/d2306a8f-59d1-4cb0-b527-ba44ca8eec3a",
        #     "http://data.nycprepared.org/ar/dataset/dycd-after-school-programs-family-support-programs-for-seniors/resource/493f52a4-0a49-4f5f-8937-78e69fb77852",
        #     "https://nycdatastables.s3.amazonaws.com/2013-07-02T15:29:20.692Z/agency-service-center-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-13T18:39:44.536Z/nyc-2012-farmers-market-list-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-10-18T21:14:52.348Z/nyc-grocery-stores-final.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-18T14:29:37.626Z/subway-entrances-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-04T17:58:59.335Z/map-of-monuments-results.csv",
        #     "https://nycdatastables.s3.amazonaws.com/2013-06-18T20:17:34.010Z/nyc-landmarks-results.csv"
            ]

    POIs = []

    for url in urls:
        csv = readCSV(url)
        print url
        for line in csv:
            latitude = line.get('latitude', None)
            longitude = line.get('longitude', None)
            if latitude is not None and longitude is not None:
                POIs.append((float(longitude), float(latitude)))
                # POIs.append((float(latitude), float(longitude)))
    return POIs


def POIsInRegion(regions, POIs):
    """
    Returns a dictionary formed by the id of a region and the number of POIs that falls in
    this region.
    """
    dict = {}

    for key, value in regions.iteritems():
        dict[key] = 0
        polygon = shape(value)
        for p in POIs:
            point = Point(p[0], p[1])
            # print point.wkt
            if polygon.contains(point):
                dict[key] += 1

    return dict


if __name__ == '__main__':
    # Geographical Features
    regions_bbox = getRegions()
    regions_number = len(regions_bbox)
    print "Regions: ", regions_number

    print "Reading POIs..."
    POIs = getPOIs()
    print len(POIs)
    print "Done Reading POIs"

    print "Calculating POIs per Region"
    POIsPerRegion = POIsInRegion(regions_bbox, POIs)
    for k,v in POIsPerRegion.iteritems():
        print k,v

Reversed the order in getPOIs function and tweaked the dict at the end.

Here's my result

{0: 9, 1: 0, 2: 2, 3: 0, 4: 4, 5: 0, 6: 0, 7: 0, 8: 13, 9: 0, 10: 3, 11: 0, 12: 0, 13: 0, 14: 1, 15: 0, 16: 0, 17: 0, 18: 0, 19: 0, 20: 0, 21: 0, 22: 0, 23: 0, 24: 0, 25: 0, 26: 0, 27: 0, 28: 0, 29: 0, 30: 0, 31: 0, 32: 0, 33: 0, 34: 0, 35: 0, 36: 0, 37: 0, 38: 0, 39: 0, 40: 0, 41: 0, 42: 2, 43: 0, 44: 0, 45: 0, 46: 0, 47: 0, 48: 0, 49: 0, 50: 0, 51: 0, 52: 0, 53: 0, 54: 0, 55: 0, 56: 0, 57: 0, 58: 0, 59: 0, 60: 0, 61: 0, 62: 0, 63: 0, 64: 0, 65: 0, 66: 0, 67: 0, 68: 0, 69: 0, 70: 0, 71: 0, 72: 0, 73: 0, 74: 0, 75: 6, 76: 0, 77: 0, 78: 1, 79: 0, 80: 0, 81: 0, 82: 0, 83: 0, 84: 2, 85: 0, 86: 0, 87: 0, 88: 0, 89: 0, 90: 0, 91: 0, 92: 0, 93: 1, 94: 0, 95: 1, 96: 1, 97: 11, 98: 1, 99: 0, 100: 1, 101: 0, 102: 3, 103: 0, 104: 3, 105: 21, 106: 0, 107: 2, 108: 1, 109: 0, 110: 3, 111: 1, 112: 0, 113: 0, 114: 0, 115: 0, 116: 0, 117: 0, 118: 1, 119: 0, 120: 0, 121: 6, 122: 0, 123: 0, 124: 2, 125: 0, 126: 2, 127: 1, 128: 2, 129: 1, 130: 0, 131: 0, 132: 0, 133: 1, 134: 8, 135: 1, 136: 0, 137: 2, 138: 1, 139: 4, 140: 3, 141: 0, 142: 0, 143: 0, 144: 0, 145: 0, 146: 0, 147: 0, 148: 0}

Best Answer

Related Solutions

Polygon – Best Way to Determine if a Point is Within a Certain Distance of a GEOJSON Polygon

[GIS] Point in Polygon – geojson – using Shapely / Python returning incorrect results

Related Question