Using Python I'm attempting to determine the state a coordinate pair resides within.

I've downloaded the US Census State And Equivalent dataset, extracted the shapefiles, and put them into Shapely objects.

sf = shapefile.Reader('shapes/tl_2017_us_state.shp')

for record in sf.shapeRecords():
    shapely_shape = shape(record.shape)

    # Ive also tried to parse the shapefile like this according to:
    # https://gis.stackexchange.com/questions/113799/how-to-read-a-shapefile-in-python?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa
    # shapely_shape = shape(record.shape.__geo_interface__)

    if isinstance(shapely_shape, Polygon):
        shapely_shape = MultiPolygon([shapely_shape])

    point = Point(37.907478, -122.369055)
    if shapely_shape.touches(point):

According to Google Maps the point (37.907478, -122.369055) should be on the California coast, however the .touches() doesn't appear to show the intersection. My suspicion is that I'm losing some accuracy when going from a shapefile to a shapely object, so I'm looking for some insight into what I'm doing wrong.

  • 1
    more likely a mismatch between where Google got it's coastline from and the Census files. Try using a buffer or dwithin rather than the very exact touches.
    – Ian Turton
    Apr 24 '18 at 12:32

You should use (longitude, latitude) syntax. point = Point(-122.369055, 37.907478)

touches returns True if the objects have at least one point in common and their interiors do not intersect with any part of the other. Overlapping features do not therefore touch. Reference: Binary Predicates

Point(-122.369055, 37.907478) is within California. Depending on your data, you should consider using overlap/within/intersects instead of touches.

enter image description here

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.