Check if a geopoint with latitude and longitude is within a shapefile

Another option is to use Shapely (a Python library based on GEOS, the engine for PostGIS) and Fiona (which is basically for reading/writing files):

import fionaimport shapelywith fiona.open("path/to/shapefile.shp") as fiona_collection:    # In this case, we'll assume the shapefile only has one record/layer (e.g., the shapefile    # is just for the borders of a single country, etc.).    shapefile_record = fiona_collection.next()    # Use Shapely to create the polygon    shape = shapely.geometry.asShape( shapefile_record['geometry'] )    point = shapely.geometry.Point(32.398516, -39.754028) # longitude, latitude    # Alternative: if point.within(shape)    if shape.contains(point):        print "Found shape for point."

Note that doing point-in-polygon tests can be expensive if the polygon is large/complicated (e.g., shapefiles for some countries with extremely irregular coastlines). In some cases it can help to use bounding boxes to quickly rule things out before doing the more intensive test:

minx, miny, maxx, maxy = shape.boundsbounding_box = shapely.geometry.box(minx, miny, maxx, maxy)if bounding_box.contains(point):    ...

Lastly, keep in mind that it takes some time to load and parse large/irregular shapefiles (unfortunately, those types of polygons are often expensive to keep in memory, too).

This is an adaptation of yosukesabai's answer.

I wanted to ensure that the point I was searching for was in the same projection system as the shapefile, so I've added code for that.

I couldn't understand why he was doing a contains test on ply = feat_in.GetGeometryRef() (in my testing things seemed to work just as well without it), so I removed that.

I've also improved the commenting to better explain what's going on (as I understand it).

#!/usr/bin/pythonimport ogrfrom IPython import embedimport sysdrv = ogr.GetDriverByName('ESRI Shapefile') #We will load a shape fileds_in = drv.Open("MN.shp")    #Get the contents of the shape filelyr_in = ds_in.GetLayer(0)    #Get the shape file's first layer#Put the title of the field you are interested in hereidx_reg = lyr_in.GetLayerDefn().GetFieldIndex("P_Loc_Nm")#If the latitude/longitude we're going to use is not in the projection#of the shapefile, then we will get erroneous results.#The following assumes that the latitude longitude is in WGS84#This is identified by the number "4326", as in "EPSG:4326"#We will create a transformation between this and the shapefile's#project, whatever it may begeo_ref = lyr_in.GetSpatialRef()point_ref=ogr.osr.SpatialReference()point_ref.ImportFromEPSG(4326)ctran=ogr.osr.CoordinateTransformation(point_ref,geo_ref)def check(lon, lat):    #Transform incoming longitude/latitude to the shapefile's projection    [lon,lat,z]=ctran.TransformPoint(lon,lat)    #Create a point    pt = ogr.Geometry(ogr.wkbPoint)    pt.SetPoint_2D(0, lon, lat)    #Set up a spatial filter such that the only features we see when we    #loop through "lyr_in" are those which overlap the point defined above    lyr_in.SetSpatialFilter(pt)    #Loop through the overlapped features and display the field of interest    for feat_in in lyr_in:        print lon, lat, feat_in.GetFieldAsString(idx_reg)#Take command-line input and do all thischeck(float(sys.argv[1]),float(sys.argv[2]))#check(-95,47)

This site, this site, and this site were helpful regarding the projection check. EPSG:4326

Here is a simple solution based on pyshp and shapely.

Let's assume that your shapefile only contains one polygon (but you can easily adapt for multiple polygons):

import shapefilefrom shapely.geometry import shape, Point# read your shapefiler = shapefile.Reader("your_shapefile.shp")# get the shapesshapes = r.shapes()# build a shapely polygon from your shapepolygon = shape(shapes[0])    def check(lon, lat):    # build a shapely point from your geopoint    point = Point(lon, lat)    # the contains function does exactly what you want    return polygon.contains(point)