point_to_raster.py

'''
Script to map points to administrative units defined by shapefiles

Requires shapefiles of administrative units from GADM (http://gadm.org/).

Alex Rutherford 2016
'''

import scipy.spatial
import shapefile,pickle,random,re,collections,csv
import numpy as np
import matplotlib.patches
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import matplotlib.path as mplPath
import time

###################
def getMid(bbox):
###################
    '''
    Helper function to get midpoint of a bounding box
    '''
    x=bbox[0]+((bbox[2]-bbox[0])/2)
    y=bbox[1]+((bbox[3]-bbox[1])/2)
    return (x,y)
###################
def getLatLong(xy):
###################
    '''
    Helper function to convert between raster indicies to lat long
    '''
    x=xy[1]
    y=xy[0]
    return (x0+(x*xSize),y0+((nRows-y)*ySize))
###################
def testPoint(point,guess,shapes0,paths0,shapes2,paths2,verbose=False):
###################
    '''
    Takes (x,y) tuple and returns index of voronoi, if hit is detected
    '''

    inCountry=False

    # First test to see if point lies within country
    for n in range(len(shapes0)):
        for p in paths0[n]:
            if p.contains_point(point):
                inCountry=True

    if not inCountry:
        if verbose:print '\tFound outside country'
        return None,None

    # If a guess of admin is provided, next test this
    if guess:
        for p in paths2[guess]:
            if p.contains_point(point):
                return guess,guess

    # If not found in test, sort admins by distance to centroid
    # and test in order

    centroidSort=collections.Counter()

    for n in xrange(len(shapes2)):
        centroidSort[n]=scipy.spatial.distance.euclidean(point,getMid(shapes2[n].bbox))

    for n,nshp in enumerate([k for k,v in centroidSort.most_common()[::-1]]):
        for p in paths2[nshp]:
            if p.contains_point(point):
                if verbose:print '\tFound in %dth shape' % n
                return nshp,nshp

    return None,None

##############
def main():
##############

    # Read in coords for country and admin 2
    # We store shapes<n> and paths<n>
    # for all admin level n

    sf = shapefile.Reader("BRA_adm2.shp")
    shapes2=sf.shapes()

    sf = shapefile.Reader("BRA_adm0.shp")
    shapes0=sf.shapes()

    paths2=collections.defaultdict(list)

    for nshp in xrange(len(shapes2)):
        pts     = np.array(shapes2[nshp].points)
        prt     = shapes2[nshp].parts
        par     = list(prt) + [pts.shape[0]]
        for pij in xrange(len(prt)):
            paths2[nshp].append(mplPath.Path(pts[par[pij]:par[pij+1]]))

    print 'Got %d paths' % len(paths2.items())

    paths0=collections.defaultdict(list)

    for nshp in xrange(len(shapes0)):
        pts     = np.array(shapes0[nshp].points)
        prt     = shapes0[nshp].parts
        par     = list(prt) + [pts.shape[0]]
        for pij in xrange(len(prt)):
            paths0[nshp].append(mplPath.Path(pts[par[pij]:par[pij+1]]))

    # Test that midpoint falls within country
    guess=None

    n,guess=testPoint((getMid(shapes0[0].bbox)),guess,shapes0,paths0,shapes2,paths2)
    print n,guess

    bbox=shapes0[0].bbox

    print bbox

    start=time.time()

    # Test 100 random points from within the bounding box

    for x,y in zip(np.random.uniform(bbox[0],bbox[1],100),np.random.uniform(bbox[2],bbox[3],100)):
        n,guess=testPoint((x,y),guess,shapes0,paths0,shapes2,paths2)

        print 'Testing (%.2f,%.2f) ' % (x,y)
        if n:
            print ': %d' % (n)
        else:
            print 'Not Found'

    end=time.time()

    print 'This took %d' % (end-start)

if __name__=='__main__':
    main()