Shapely not able to precisely find points inside polygon

Question

I have several points and I would like to check if they are contained on a polygon. The polygon and the points are represented by latitudes and longitudes.

Following is the code of to reproduce my scenario and the Google Maps print screen of what it looks like the polygon, the points inside/outside the polygon as per Shapely.

import pyproj
from shapely.geometry import Polygon, Point
from shapely.ops import transform
from functools import partial
import numpy as np

polygon_array = [(1.4666748046875, 49.088257784724675),
                 (1.4447021484375, 47.42808726171425),
                 (2.889404296875, 47.42808726171425),
                 (2.8729248046875, 49.08466020484928),
                 (-0.0054931640625, 47.97521412341619),
                 (0.010986328125, 46.18743678432541),
                 (1.4227294921875, 46.1912395780416),
                 (1.4337158203125, 48.02299832104887),
                 (-1.043701171875, 46.65320687122665),
                 (-1.043701171875, 44.6061127451739),
                 (0.0164794921875, 44.5982904898401),
                 (-0.0054931640625, 46.6795944656402)]

simple_polygon = Polygon(polygon_array)
projection = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'),
    pyproj.Proj('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs'))
polygon = transform(projection, simple_polygon)

for latitude in np.arange(44.5435052132, 49.131408414, 0.071388739257):
    for longitude in np.arange(-0.999755859375, 2.99926757812, 0.071388739257):
        point = transform(projection, Point(longitude, latitude))
        if polygon.contains(point):
            print "%s, %s" % (latitude, longitude)

Here is what the polygon looks on a map:

Here is what looks like the points (here represented as markers) "inside" the polygon:

And the points "outside":

The problem here is that the points are clearly way off the polygon, inside or outside. I am new to this projection scheme so I may be missing something.

Thank you in advance

Answer 1

Your polygon doesn't look anything like the picture you drew (best I can tell).

fiddle

code snippet:

function initialize() {
  var map = new google.maps.Map(
    document.getElementById("map_canvas"), {
      center: new google.maps.LatLng(37.4419, -122.1419),
      zoom: 13,
      mapTypeId: google.maps.MapTypeId.ROADMAP
    });
  var polygon_array = [{
    lng: 1.4666748046875,
    lat: 49.088257784724675
  }, {
    lng: 1.4447021484375,
    lat: 47.42808726171425
  }, {
    lng: 2.889404296875,
    lat: 47.42808726171425
  }, {
    lng: 2.8729248046875,
    lat: 49.08466020484928
  }, {
    lng: -0.0054931640625,
    lat: 47.97521412341619
  }, {
    lng: 0.010986328125,
    lat: 46.18743678432541
  }, {
    lng: 1.4227294921875,
    lat: 46.1912395780416
  }, {
    lng: 1.4337158203125,
    lat: 48.02299832104887
  }, {
    lng: -1.043701171875,
    lat: 46.65320687122665
  }, {
    lng: -1.043701171875,
    lat: 44.6061127451739
  }, {
    lng: 0.0164794921875,
    lat: 44.5982904898401
  }, {
    lng: -0.0054931640625,
    lat: 46.6795944656402
  }];
  for (var i = 0; i < polygon_array.length; i++) {
    var marker = new google.maps.Marker({
      map: map,
      position: polygon_array[i],
      title: "" + i
    })
  }
  var polygon = new google.maps.Polygon({
    map: map,
    paths: [polygon_array],
    fillOpacity: 0.5,
    strokeWeight: 2,
    strokeOpacity: 1,
    strokeColor: "red",
    fillColor: "red"
  });
  var bounds = new google.maps.LatLngBounds();
  for (var i = 0; i < polygon.getPaths().getAt(0).getLength(); i++) {
    bounds.extend(polygon.getPaths().getAt(0).getAt(i));
  }
  map.fitBounds(bounds);
}
google.maps.event.addDomListener(window, "load", initialize);

html,
body,
#map_canvas {
  height: 100%;
  width: 100%;
  margin: 0px;
  padding: 0px
}

<script src="https://maps.googleapis.com/maps/api/js"></script>
<div id="map_canvas"></div>