Using pysal to map location of coordinates?

Question

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)

Answer 1

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.