Read a shapefile of polygons and a list of points, and return for each point the polygon that contains it. Attemp 1, works but takes almost 1/2 second to calculate each point. Next step: Find a library to do this :-)

points_in_polygons1.py

import shapefile
import pandas as pd


# Point in polygon
def pip(point,poly):
    x = point[0]
    y = point[1]
    n = len(poly)
    inside = False
    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y
    return inside


def overlay_pip(point, poly_shapes_reader):
    for sr in poly_shapes_reader.shapeRecords():
        # if not in bounding, do not even try
        if sr.shape.bbox[0] <= point[0] <= sr.shape.bbox[2] and  sr.shape.bbox[1] <= point[1] <= sr.shape.bbox[3]:
            if pip(point, sr.shape.points):
                return sr.record
    return [0, 0, "Not found"]




if __name__ == "__main__":
    #provinces = shapefile.Reader("/data/geobrowsing/spain_gis/SHP_WGS84_PROVINCIAS/PROVINCIAS_WGS84.shp")
    provinces = shapefile.Reader("/data/geobrowsing/spain_gis/SHP_WGS84_PROVINCIAS_SIMPLE/PROVINCIAS_SIMPLE_OK.shp")

    antennas = pd.read_csv('/data/geobrowsing/ES_UCR20130926_head.csv', header=None,
                            names=['cellid', 'lat', 'lon', 'height', 'centlat', 'centlon', 'zipcode', 'tech', 'type', 'freq', 'beam', 'azimuth'])


    for x in antennas[['cellid', 'lon', 'lat']].values:
        print str(list(x[0:3]) + overlay_pip(x[1:3], provinces))

Reading from csv without pandas

import shapefile
import csv
import sys

def read_csv_manual(filename):
   fp = open(filename, 'rb')
   reader = csv.reader(fp)
   output = list()
   for row in reader:
       output.append([row[0], float(row[2]), float(row[1])])
   return output

# Point in polygon
def pip(point,poly):
    x = point[0]
    y = point[1]
    n = len(poly)
    inside = False
    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y
    return inside


def overlay_pip(point, poly_shapes_reader):
    for sr in poly_shapes_reader.shapeRecords():
        # if not in bounding, do not even try
        if sr.shape.bbox[0] <= point[0] <= sr.shape.bbox[2] and  sr.shape.bbox[1] <= point[1] <= sr.shape.bbox[3]:
            if pip(point, sr.shape.points):
                return sr.record
    return [0, 0, "Not found"]


if __name__ == "__main__":
    print sys.argv[1]
    provinces = shapefile.Reader("/home/acg/geobrowsing/SHP_WGS84_PROVINCIAS_SIMPLE/PROVINCIAS_SIMPLE_OK.shp")

    ants = read_csv_manual('/home/acg/geobrowsing/UCR/x' + sys.argv[1])

    for x in ants:
        print str(list(x[0:3]) + overlay_pip(x[1:3], provinces))