0

I have a list of points describing the boundaries of Spain. I want to be able to tell whether a pair of lat,lon is within these boundaries. I have tried the following:

import shapefile
import matplotlib.pyplot as plt
from shapely.geometry import MultiPoint, Point, Polygon
from shapely.geometry.polygon import Polygon
sf = shapefile.Reader(r"\ESP_adm0.shp")
shapes = sf.shapes()
lat = []; lon = []
for i in range(len(shapes[0].points)):
    lon.append(shapes[0].points[i][0]);lat.append(shapes[0].points[i][1])

I know I am retrieving the points, because I'm able to plot and get the desired results:

plt.plot(lon,lat,'.', ms=0.1)

(plot in the link below) plot result

I do the following to get the poitns into a polygon:

    coords = list(zip(lat,lon))
    spain_pol = Polygon(coords)

And then I use the contains function.

    spain_pol.contains(Point(0,42))
    spain_pol.contains(Point(42,0))

These both return false. In fact I haven't been able to get a single point I've tried to return a True.

I have tried all sorts of things, and I think I must be missing something fundamental.

Perhaps the facts Spain has islands and there's more than one polygon is the problem?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Alpha
  • 21
  • 1
  • 7
  • 2
    What is the projection system of ESP_adm0.shp? – neogeomat Mar 10 '18 at 11:30
  • I got the file from http://www.gadm.org/, not sure what the projection is. That said, I can see the coordinates of the points, and from the plot, I can tell that at 0,42, it should return TRUE, regardless of the projection, because in my mind the projection is prior to the coordinates, and in the plot you can see the coordinates post-projection. I could be completely wrong here, I'm a total non expert. Thank you for your help – Alpha Mar 10 '18 at 12:07

1 Answers1

2

Just in case someone else has the same issue. The following code worked perfectly:

import fiona
from shapely.geometry import MultiPoint, Point, Polygon,shape
from shapely.geometry.polygon import Polygon

multipol = fiona.open(r"C:\Users\Jordi\Downloads\ESP_adm_shp\ESP_adm0.shp")
multi = next(iter(multipol))

point = Point(0,42)
point.within(shape(multi['geometry']))

This returns a very welcome "True" :)

Alpha
  • 21
  • 1
  • 7