points_in_polygons_or_near.py

points_in_polygons_or_near.py

import shapefile
import pandas as pd
from shapely.geometry import Polygon, Point

provinces = shapefile.Reader("/data/geobrowsing/spain_gis/SHP_WGS84_PROVINCIAS_SIMPLE/PROVINCIAS_SIMPLE_OK.shp")
provs = list()

for sr in provinces.shapeRecords():
    provs.append(sr.record + [Polygon(sr.shape.points),] )
    #print sr.record
    #[2, '02', 'Albacete']
    #res = pip(test_point, sr.shape.points)
    #print str(res) + "->" + str(sr.record) + str(len(sr.shape.points)) + "__" + str(sr.shape.bbox)


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

ants = list()
for u in ucr[['cellid', 'lon', 'lat']].values:
    ants.append(list(u) + [Point(u[1:3]), ])


for a in ants:
    found = False
    for prov in provs:
        if (not found) and (a[3].within(prov[3])):
            print a[0:3] + prov[0:3] + ['inside', ]
            found = True
    if not found:  # If not found, return closest
        distances = list()
        for prov in provs:
            distances.append(a[3].distance(prov[3]))
        posmin = distances.index(min(distances))
        print a[0:3] + provs[posmin][0:3] + ['not found', ]