7

I'm trying to map coordinate data to a municipality. Reverse geocoders aren't as accurate as I need them to be (I was using geocoder pkg before), so I found a .shp that maps and names all the municipalities I need to refer to.

I have 9000 records in a CSV that need to be mapped, and was wondering If anyone knew how to check if the individual coordinates are in these shapes.

I'm not familiar with .shp files, and don't know if this would be possible through a package like pyshp


With help, I got it. Here is the code in case anyone has a similar issue in the future, it's sloppy and I may update when I'm done cleaning it up. Loading up the library of shapes is pretty slow, just to warn. And checking each point against each shape has high complexity, but for my sample of 9000 values against ~45 shapes it ran pretty quickly:

import shapefile as shp
from shapely.geometry import Point
import csv
from shapely.geometry.polygon import Polygon

fileName = 'TestLatLong'                #Open CSV file with points to check   
sf = shp.Reader("NewMuni/NewMuni111")   #Open Shapefile with shapes to check points against
sfRec = sf.records() #Read records in shapefile

n = 0
m = 1
coor = ''
coorDict = {}
matplotDict = []
muniFinal = {}

for shape in sf.shapeRecords(): #Iterate through shapes in shapefile
    x = [i[0] for i in shape.shape.points[:]] #Initially for use in matplotlib to check shapefile
    y = [i[1] for i in shape.shape.points[:]] #Initially for use in matplotlib to check shapefile
    for i in x:
        matplotDict.append((x[x.index(i)],y[x.index(i)])) #Convert coordinates to be read by Shapely pkg

    munishp = Polygon(matplotDict)
    muniFinal[sfRec[n][1]] = munishp #Store shape in dictionary with key of municipality

    matplotDict = [] #refresh coordinate store for next shape   
    n += 1 


n = 0    



with open(fileName + '.csv') as csvfile:
    reader = csv.DictReader(csvfile)
    for row in reader:

        coor = (row['latitude'],row['longitude'])
        rlat = float(row['latitude'])
        rlong = float(row['longitude'])
        if coor == ' ,   ' or coor == ', ':
            coorDict[row['PrimaryKey']] = 'No Data' #PrimaryKey is my primary key that I will use to write the data back into the .csv 
        else:
            if float(row['longitude']) > 0:
                coorDict[row['PrimaryKey']] = (rlat,rlong)
            else:
                coorDict[row['PrimaryKey']] = (rlong,rlat)
        m += 1

#proof of concept- save the results however you'd like
for j in coorDict:
    for k in muniFinal:
        if muniFinal[k].contains(Point(coorDict[j])):
            print(j, 'in', k)
3
  • Did I understand correctly that you have one .csv file with X, Y columns and one polygon shapefile with the municipalities boundaries and you would like to create a new point shapefile that would contain all the point features from .csv file with extra column MunicipalityName coming from the boundaries shapefile? Commented Jul 27, 2017 at 18:15
  • Try converting to points and running an intersect operation.
    – Aaron
    Commented Jul 27, 2017 at 18:33
  • Hi Alex, sorry for such a delayed response, but correct: I found a resource that has a folder that can be uploaded to arcGIS (.shp,shx,sbn,dbf, etc...) then I have a csv file with latitude and longitude columns for every entry. I'm trying to explore the docs on shapefile, shapely, and then use pandas and csv to append data to each row regarding which municipality (polygon) they are inside Commented Aug 1, 2017 at 14:33

2 Answers 2

8

PyShp will let you read the shapefile, but won't help you figure out if a point is in a boundary. You'll have to combine it with something like Shapely to do the geometric calculations. Luckily, the two modules can interoperate through the Python Geo Interface. Some basic functionality (untested) would be like:

import shapefile
from Shapely.geometry import Point # Point class
from Shapely.geometry import shape # shape() is a function to convert geo objects through the interface

point_to_check = (1234,5678) # an x,y tuple
shp = shapefile.Reader('path/to/shp') #open the shapefile
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()     
for i in len(all_shapes):
    boundary = all_shapes[i] # get a boundary polygon
    if Point(point_to_check).within(shape(boundary)): # make a point and see if it's in the polygon
       name = all_records[i][2] # get the second field of the corresponding record
       print "The point is in", name     

Make sure your points and boundaries are in the same coordinate system!

3
  • I had trouble getting Shapely's shape function to work, so I just used Polygon instead, but this really helped! Commented Aug 4, 2017 at 20:48
  • 1
    for those cutting-and-pasting, the actual modules are: import shapefile from shapely.geometry import Point from shapely.geometry import shape Commented Mar 4, 2018 at 22:21
  • Just one general point for working on such problems. Point gets (longitude, latitude), not vice versa.
    – mah65
    Commented Nov 8, 2020 at 0:50
0

Following Marc's answer, it should be

for i in range(1,len(all_shapes)):
    boundary = all_shapes[i] # get a boundary polygon
    #print(shape(boundary))
    if Point(point_to_check).within(shape(boundary)): # make a point and see if it's in the polygon
       c = all_records[i][2] # get the second field of the corresponding record
       print(all_records[i][2])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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