def pointsToVoronoiGridShapefile(lat, lon, vor_shp_path, extent=None):
    """
    Converts points to shapefile grid via voronoi
    """
    voronoi_centroids = _get_voronoi_centroid_array(lat, lon, extent)

    # set-up output polygon shp
    log("Creating output polygon shp {0}"
        .format(os.path.basename(vor_shp_path)))
    if os.path.exists(vor_shp_path):
        os.remove(vor_shp_path)
    drv = ogr.GetDriverByName('ESRI Shapefile')
    outShp = drv.CreateDataSource(vor_shp_path)
    osr_geographic_proj = osr.SpatialReference()
    osr_geographic_proj.ImportFromEPSG(4326)
    layer = outShp.CreateLayer('', osr_geographic_proj, ogr.wkbPolygon)
    layer.CreateField(ogr.FieldDefn('GRID_LAT', ogr.OFTReal))
    layer.CreateField(ogr.FieldDefn('GRID_LON', ogr.OFTReal))
    layerDefn = layer.GetLayerDefn()

    # find nodes surrounding polygon centroid
    # sort nodes in counterclockwise order
    # create polygon perimeter through nodes
    log("Building Voronoi polygons...")
    # compute voronoi
    voronoi_manager = Voronoi(voronoi_centroids)
    voronoi_vertices = voronoi_manager.vertices
    voronoi_regions = voronoi_manager.regions
    for point_id, region_index in enumerate(voronoi_manager.point_region):
        vert_index_list = np.array(voronoi_regions[region_index])
        voronoi_centroid = voronoi_centroids[point_id]
        voronoi_poly_points = _get_voronoi_poly_points(vert_index_list,
                                                       voronoi_vertices,
                                                       voronoi_centroid)
        if len(voronoi_poly_points) == 4:
            poly = ogr.Geometry(ogr.wkbPolygon)
            ring = ogr.Geometry(ogr.wkbLinearRing)
            for node in voronoi_poly_points:
                ring.AddPoint(node[0], node[1])

            # grab first node to close ring
            ring.AddPoint(voronoi_poly_points[0][0], voronoi_poly_points[0][1])

            poly.AddGeometry(ring)
            feat = ogr.Feature(layerDefn)
            feat.SetField('GRID_LON', float(voronoi_centroid[0]))
            feat.SetField('GRID_LAT', float(voronoi_centroid[1]))
            feat.SetGeometry(poly)
            layer.CreateFeature(feat)
Converts points to shapefile grid via voronoi