We help IT Professionals succeed at work.

# How do I search a shapefile using GPS coordinates?

on
Medium Priority
132 Views
I have a standard shapefile.  It contains many shapes with their attributes.
I have a list of GPS coordinates that I know for sure are within the area of that shapefile.
How -- preferably with Python -- can I supply the lat/long and get an handle to that shape?
Comment
Watch Question

## View Solution Only

HelpDesk / Programmer
CERTIFIED EXPERT

Commented:
Hi,

you could try this python libraries:
1.  pyshp
2.  pyshapelib
3.  shapelib

Hope it helps,
Istvan

Commented:
Thank you, I've tried 1 but cannot see how I could find the index of a shape (record) by supplying a GPS point (lat/long).
HelpDesk / Programmer
CERTIFIED EXPERT

Commented:
i never used this library, but i made a simple app in python, hope it gives you some idea:

code:
``````import shapefile

# print(shp)

shapes = sf.shapes()
pointindex=-1
shapeIndex=-1
for shape in shapes:
shapeIndex=shapeIndex+1
for point in shape.points:
pointindex=pointindex+1
#print('shape[idx:',shapeIndex,']= X(',point,'),Y(',point,')','\n')
if (point==-122.465628 and point==37.66665) or (point==-122.469816 and point==37.675509):
print('Found Coordinate ! Shape index: ', shapeIndex,', Point index: ',pointindex)
``````
------------------------------------------------------
you can try it in the online compiler:

hope it helps
Istvan
HelpDesk / Programmer
CERTIFIED EXPERT

Commented:
the shape file in google maps looks like this:

Commented:
Excellent!  The code gives me all the points of the shape.  Providing the lat/long of any of those points identifies the shape.
The question is how do I know if the lat/long is within the polygon defined by the lines connecting those points?
HelpDesk / Programmer
CERTIFIED EXPERT
Commented:
From what i read: pyshp library can read the shape file, but can't figure out if a point is in a boundary.
To make this kind of calculation we need use Shapely. Combining the two we can read and calculate.

The points and boundaries should be in same coordinate system.

you could try like this:
``````import shapefile
from shapely.geometry import Point
from shapely.geometry import shape

point = (-122.44553987365936,37.68154318087043) # an x,y tuple
shp = shapefile.Reader('shapefiles/blockgroups.shp') #open the shapefile
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()
for i in range(len(all_shapes)):
boundary = all_shapes[i] # get a boundary polygon
if Point(point).within(shape(boundary)): # make a point and see if it's in the polygon
pop1990 = all_records[i] # get the second field of the corresponding record
print("---------------------------------------------------------------")
print("Point(",point,",",point,") is in boundary: ",Point(point).within(shape(boundary))," / Shape index: ",i,"  :")
print("Secound field value: ",pop1990)
print("---------------------------------------------------------------")
for j in range(len(all_records[i])):
print("Fields[",j,"]:",all_records[i][j]);
``````

try the code: online code

Similar question:
1. Link -( the code above is based on this info )