Fiona and Shapely are better suited to this task, which is essentially a spatial join (contains). Here is some code that does what you want. However, I will caution you that this is not the best way to do this (algorithmically), and if I were you, I would be using a spatial database to complete the join. I'll demonstrate that, too.
import fiona # Reads and writes spatial data files
from shapely.geometry import Point, shape # Does geometry!
# shapefile
COUNTY_FILE = './English_FILES/English Ceremonial Counties.shp'
# Open text file of coordinates
with open('./sundayCoordinates.txt') as coords:
# Clean the text points
coords = [g.replace('(','').replace(')','').strip().split(',') for g in coords]
points = [Point(float(g[1]), float(g[0])) for g in coords] # shapely Points
with fiona.open(COUNTY_FILE, 'r') as counties:
for county in counties:
g = county['geometry'] # Geometry of a county
try:
assert g['type'] == "Polygon" # Sanity check that it is indeed a Polygon
except AssertionError:
print "You have a problem with %s's geometry" % county['properties']['NAME']
continue # Move onto next county as this one is invalid
g = shape(g) # shapely shape from geometry
for p in points: # step through coordinates (shapely Points)
if g.contains(p):
print "%s is in %s" % (p, county['properties']['NAME'])
# break # uncomment this command if you just want to illustrate output
else:
continue # Move on to next point
Uncommenting the break command to prevent full execution, you get the following output:
POINT (-0.42417508 54.00307178) is in East Riding of Yorkshire
POINT (-0.63751549 52.65048587) is in Rutland
POINT (-2.73563043 52.7027027) is in Shropshire
POINT (-1.1853677 51.4669939) is in Berkshire
POINT (-1.8267535 51.6049064) is in Wiltshire
You have a problem with Isle of Wight's geometry
POINT (-1.73318956 52.45102757) is in West Midlands
POINT (-1.86776924 53.73871304) is in West Yorkshire
POINT (-3.0189 53.4549) is in Merseyside
POINT (-2.2516401 53.4868507) is in Greater Manchester
POINT (-1.38795678 54.90352965) is in Tyne and Wear
POINT (-1.4647484 53.3765869) is in South Yorkshire
POINT (-0.533208 51.9470787) is in Bedfordshire
POINT (-2.43301135 53.08910114) is in Cheshire
POINT (-2.6281595 51.4289207) is in County of Bristol
POINT (-2.80496618 52.158666) is in Herefordshire
POINT (-2.13333 52.9) is in Staffordshire
POINT (0.91444444 52.33555556) is in Suffolk
POINT (-0.412932 51.403928) is in Surrey
POINT (-1.26114 52.3675) is in Warwickshire
POINT (-0.319 51.066) is in West Sussex
POINT (-2.502849 52.3154415) is in Worcestershire
POINT (0.00315428 51.50285246) is in Greater London
You have a problem with Devon's geometry
You have a problem with Northumberland's geometry
You have a problem with Cornwall's geometry
POINT (-0.56444444 51.59861111) is in Buckinghamshire
POINT (0.138 52.194) is in Cambridgeshire
POINT (-2.59899508 54.64857178) is in Cumbria
POINT (-2.008166 53.371521) is in Derbyshire
POINT (-2.4697236 50.6370497) is in Dorset
POINT (0.37338294 50.81937006) is in East Sussex
POINT (0.73111111 51.86805556) is in Essex
POINT (-2.09135056 51.89923046) is in Gloucestershire
You have a problem with Hampshire's geometry
POINT (-0.39032 51.656489) is in Hertfordshire
POINT (0.27227 51.43914) is in Kent
POINT (-1.17675776 52.45203366) is in Leicestershire
POINT (-0.53377 53.22462) is in Lincolnshire
POINT (1.29620869 52.62843255) is in Norfolk
POINT (-0.8871239 52.051104) is in Northamptonshire
POINT (-1.0783 53.9454193) is in North Yorkshire
POINT (-0.95 53.0667) is in Nottinghamshire
POINT (-1.2328819 51.7248487) is in Oxfordshire
You have a problem with Somerset's geometry
POINT (-2.68847466 53.772178) is in Lancashire
POINT (-1.3139652 54.56139774) is in Durham
The problematic geometries are (I'm guessing) because they are multipolygons. I think you can probably safely use contains
with these geometries: try it out!
A better way, in my opinion—and this is totally optional—is to use PostGIS. I imported your counties shapefile to PostGIS with shp2pgsql
: shp2pgsql English_FILES/English\ Ceremonial\ Counties.shp public.counties | psql -h localhost -d my_database -U postgres
. I also did the same with your coordinates text file. Then your code to match the two datasets simply becomes:
SELECT ST_AsText(p.geom), s.name
FROM coords AS p, counties AS s
WHERE ST_Within(p.geom, s.geom)
And you get, for example:
"POINT(0.00315428 51.50285246)";"Greater London"
"POINT(-1.86776924 53.73871304)";"West Yorkshire"
"POINT(-0.21083 51.4955)";"Greater London"
"POINT(-0.15269279 51.51435128)";"Greater London"
"POINT(-0.16515719 51.49203227)";"Greater London"
"POINT(-0.13005146 51.59350874)";"Greater London"
PostGIS is amazingly powerful, so I recommend that ultimately. Happy learning.
points
variable, and then never use it. I don't understand that the point of yourfor line in coordinates
iteration is. – alphabetasoup Dec 4 '15 at 10:43coordinates in counties
is False. If you make your data available somewhere I could take a more thorough look, but start by adding a bunch of print statements to better understand what is happening. For instance, add anelse
to theif
inside thewhile
loop, andprint
something. You'll see the infinite loop. – alphabetasoup Dec 4 '15 at 10:57points in counties
it doesn't seem to work, even though I know the points are in the uk? I uploaded the code I am using, the coordinates and the shapefile to dropbox, if you could take a quick look that would be really appreciated: dropbox.com/sh/52ktt4xp6goac5g/AADn6XNM_aCPzj3K0vmY7C4ba?dl=0 – JTH Dec 4 '15 at 11:09in
compares, essentially, that a chosen element is within a set. E.g.'a' in ['a','b','c'] is True
yet'd'
is not in that same set. It is not performing a geometric contains operation. – alphabetasoup Dec 4 '15 at 11:27