1

I have a shapefile that breaks the city of San Francisco up into a number of different districts depending on who your supervisor is.

I need to make a lookup tool that takes a single long lat value and says which supervisor district the point lies in. The supervisor district is one of the attributes of the shapefile.

How would I go about doing this in an efficient manner?

Edit 1: Seems like something like this is what I'm looking for. Trying to figure out how I access the attributes once an intersection is found

  • This is a simple problem (ArcGIS solution), but we need to know some more information (e.g., what software you are using). Welcome to GIS.SE! – Paul Sep 28 '15 at 22:28
  • What software or API are you using? Most GIS libraries have an intersect select by location or within function which can be used to decide what polygon the point is within. As Paul said this is really simple to lookup. – Michael Stimson Sep 28 '15 at 22:29
  • I need to do this in python with an open source solution. Trying to see right now if I can do it using shapely and fiona. – Luke Sep 28 '15 at 22:30
  • I've not had any experience with either of those but OGR certainly does, gdal.org/…. The most basic way is first create your point as an OGRgeometry then iterate each polygon in your shapefile (OGRlayer) until within is true; beware spatial references - make sure your coordinates for your point match your polygons. OGR is part of GDAL, open source, available in Windows, Linux and Mac for languages python, C, C++, C#, VB, ruby and many more. Is this server/web or to be run locally? – Michael Stimson Sep 28 '15 at 23:29
2

I ended up doing it this way. The shapefile was too small for an rtree to make any difference in speed:

import fiona
from shapely.geometry import Point, shape

def long_lat_to_district(longit, lat, shape_file):

    mypoint = Point(longit, lat)   

    with fiona.open(shape_file) as shp:
        polygons = [poly for poly in shp]

    poly_idx = [i for i, poly in enumerate(polygons)
                if mypoint.within(shape(poly['geometry']))]

    if not poly_idx:
        return None
    else:
        # Take first polygon that overlaps since may overlap with several if on border
        match = polygons[poly_idx[0]]
        return match['properties']

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.