Profile between two points on Map Python

#!/usr/bin/python
import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt
from matplotlib.mlab import griddata
from scipy.interpolate import griddata
import matplotlib.mlab as ml
import numpy
import pylab
from scipy.signal import convolve2d


#-- Generate some data...

def moving_average_2d(data, window):
    """Moving average on two-dimensional data.
    """
    # Makes sure that the window function is normalized.
    window /= window.sum()
    # Makes sure data array is a numpy array or masked array.
    if type(data).__name__ not in ['ndarray', 'MaskedArray']:
        data = numpy.asarray(data)
    # The output array has the same dimensions as the input data
    # (mode='same') and symmetrical boundary conditions are assumed
    # (boundary='symm').
    return convolve2d(data, window, mode='same', boundary='symm')

######## extrapolate NaN
def extrapolate_nans(x, y, v):
    if np.ma.is_masked(v):
        nans = v.mask
    else:
        nans = np.isnan(v)
    notnans = np.logical_not(nans)
    v[nans] = scipy.interpolate.griddata((x[notnans], y[notnans]), v[notnans],
        (x[nans], y[nans]), method='nearest').ravel()
    return v


### input data
data = np.genfromtxt('input2d.txt')
x1 = data[:,0]
y1 = data[:,1]
z1 = data[:,2]

### gridding
numcols, numrows = 50, 50
numel=numcols*numrows
xi = np.linspace(min(x1), max(x1), numcols)
yi = np.linspace(min(y1), max(y1), numrows)
xi, yi = np.meshgrid(xi, yi)
x, y, z = x1, y1, z1
zi = ml.griddata(x, y, z, xi, yi, interp='nn')
extrapolate_nans(xi,yi,zi)

### window for moving average
m, n = 5, 5    # The shape of the window array
# Calculating a couple of smoothed data.
win = numpy.ones((m, n))
#zis = ndimage.gaussian_filter(zi, sigma=25.0, order=0)
zis = moving_average_2d(zi, win)




### Make a line with "num" points...
x0, y0 =820000, 3200000
x1, y1 =900000, 3280000

###convert  coordinate into pixel
xmin=min(xi.reshape((numel,1)))
xmax=max(xi.reshape((numel,1)))
ymin=min(yi.reshape((numel,1)))
ymax=max(yi.reshape((numel,1)))

x0p=numcols*(x0-xmin)/(xmax-xmin)
x1p=numcols*(x1-xmin)/(xmax-xmin)
y0p=numrows*(y0-ymin)/(ymax-ymin)
y1p=numrows*(y1-ymin)/(ymax-ymin)
num = 100

x2p, y2p = np.linspace(x0p, x1p, num), np.linspace(y0p, y1p, num)  ###pixel coordinate
x2, y2 = np.linspace(x0, x1, num), np.linspace(y0, y1, num)        ###original coordinate

###Extract the values along the line, using cubic interpolation (pixel)
xyz = scipy.ndimage.map_coordinates(zi, np.vstack((x2p,y2p)))
xyzs = scipy.ndimage.map_coordinates(zis, np.vstack((x2p,y2p)))

###plot
fig=plt.figure()
ax=fig.add_subplot(2,1,1)
plt.tricontourf(xi.ravel(), yi.ravel(),zis.ravel(),50)
plt.plot([x0, x1], [y0, y1], 'ro-')
plt.axis((min(xi.ravel()),max(xi.ravel()),min(yi.ravel()),max(yi.ravel())))

plt.colorbar()

ax=fig.add_subplot(2,1,2)
plt.plot(x2,xyz,'r')
plt.plot(x2,xyzs,'b')
plt.show()