# imports
    from osgeo import gdal, ogr, osr
    import os
    import matplotlib.pyplot as plt
    import scipy.interpolate as il
    import numpy as np
    from itertools import cycle


    # obtaining xyz to interpolate
    a=[]        # to store all x valuesa
    b=[]        # to store all corresponding y values
    c=[]        # to store all corresponding z values
    # input output
    shapefile = r"P:May2014TestLielaisIV_Bp_Schichten.shp"
    outras=r"P:Interpolatio_spline_with_bariersoutras.tif"
    ## opening a shape file
    DriverName = ogr.GetDriverByName("ESRI Shapefile")
    dataSource = DriverName.Open(shapefile, 0)
    print dataSource
    layer = dataSource.GetLayer()
    srs= layer.GetSpatialRef()
    #looping through features to get geometry reference for each ........
    for i in range(0, layer.GetFeatureCount()):
        # Get the input Feature
        inFeature = layer.GetFeature(i)
        geom = inFeature.GetGeometryRef()
        x= geom.GetX()
        a.append(x)
        y= geom.GetY()
        b.append(y)
        z= inFeature.GetField("of_wth_01")
        c.append(z)
    # obtaining max and min extent
    xmin,xmax,ymin,ymax = [min(a),max(a),min(b),max(b)]

    #size of 1 m grid
    nx = (int(xmax - xmin)+1)
    ny = (int(ymax - ymin)+1)
    # Generate a regular grid to interpolate the data.
    xi = np.linspace(xmin, xmax, nx)#, endpoint=False)
    yi = np.linspace(ymin, ymax, ny)
    xi, yi = np.meshgrid(xi, yi)
    # Interpolate the values of z for all points in the rectangular grid
    zi = il.griddata((x, y), c, (xi, yi),method='linear')

# Interpolate the values of z for all points in the rectangular grid
zi = il.griddata((x, y), c, (xi, yi),method='linear')
#                 ^  ^