3

I am trying to map which county a set of coordinates I have are in. I want the code to go through each coordinate, check which county it is in, and then print the corresponding county. However, currently when ran the cursor in idle just flashes and nothing appears to happen. I'm not sure if the code is correctly checking the counties or whether it can get the the names of the counties from the shapefile?

import pysal as ps
from ast import literal_eval

COUNTIES = '/Users/JoshuaHawley/Desktop/counties/English Ceremonial Counties.shp'
shps = ps.open(COUNTIES)

counties = shps.read()

f = open('/Users/JoshuaHawley/sundayCoordinates.txt').read()
seq = literal_eval(f.replace("\n", ","))

coordinates = seq[0:]
for line in coordinates:
    points = coordinates[:2]

NotFound = True

while NotFound == True:
    if coordinates in counties:
        NotFound = False
        name_of_county = counties.getNamePolygon()
        print(name_of_county)
  • You set the points variable, and then never use it. I don't understand that the point of your for line in coordinates iteration is. – alphabetasoup Dec 4 '15 at 10:43
  • I wasn't sure how to iterate over one coordinate at a time? I want it to take in one coordinate and check its location, before then moving on to the next one rather than importing all the coordinates at once? Very new to python, so may be making a mistake. – JTH Dec 4 '15 at 10:50
  • Yes, there are some major syntax and logic problems. The reason you are getting no feedback is that your while loop never breaks because coordinates 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 an else to the if inside the while loop, and print something. You'll see the infinite loop. – alphabetasoup Dec 4 '15 at 10:57
  • @alpha-beta-soup I added the print command in and saw that it continually kept printing the same word. Even when using points 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:09
  • I might try tomorrow. However I think your biggest misunderstanding is that the reserved word in 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
2

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.

  • 1
    For Spatial join with Fiona and Shapely, look at More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc – gene Dec 4 '15 at 14:49
  • @alpha-beta-soup works really well, thanks for the help, have been able to write the code to count frequencies of county occurrences now too. Just wondering what you exactly mean by the contains function for those that do not work? – JTH Dec 4 '15 at 16:39
  • I only allowed contains to run on polygon shapes, not multipolygon shapes. It should work if to remove that check. – alphabetasoup Dec 4 '15 at 17:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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