Grid Rotation-Smoothing-Plot Python

#!/usr/bin/python
from __future__ import division
import numpy as np
import scipy
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
import scipy.ndimage as ndimage
import numpy
import pylab
from scipy.signal import convolve2d
from matplotlib import mlab, cm
from mpl_toolkits.mplot3d import Axes3D

pi     = scipy.pi
dot    = scipy.dot
sin    = scipy.sin
cos    = scipy.cos
ar     = scipy.array
rand   = scipy.rand
arange = scipy.arange
plot   = pylab.plot
show   = pylab.show
axis   = pylab.axis
grid   = pylab.grid
title  = pylab.title


#the function
def Rotate2D(pts,cnt,ang):
    return dot(pts-cnt,ar([[cos(ang*np.pi/180.),sin(ang*np.pi/180.)],[-sin(ang*np.pi/180.),cos(ang*np.pi/180.)]]))+cnt

#the code for test

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

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

numcols, numrows = 50, 50
numel=numcols*numrows
xi = np.linspace(min(x1)-5000, max(x1), numcols)
yi = np.linspace(min(y1)+20000, max(y1)+5000, numrows)
xi, yi = np.meshgrid(xi, yi)

pts=np.column_stack((xi.reshape((numel,1)),yi.reshape((numel,1))))

xi1=xi.reshape((numel,1))
yi1=yi.reshape((numel,1))

x0=(min(x1)+max(x1))/2.
y0=(min(y1)+max(y1))/2.

ots = Rotate2D(pts,ar([x0,y0]),38)       #rotate 38 deg: positive anti clock wise  at anchor point x0,y0
x2=ots[:,0]
y2=ots[:,1]


#### QC Rotated Grid
fig=plt.figure()
xxx=plt.plot(x1,y1, 'k. ')
plt.plot(x2,y2, 'b. ')
plt.title('Grid QC: black (data) blue (rotated grid)')
plt.xlabel('X-COORD')
plt.ylabel('Y-COORD')
plt.show()
###########

xir=x2.reshape((numcols,numrows))
yir=y2.reshape((numcols,numrows))
x, y, z = x1, y1, z1
zir = ml.griddata(x, y, z, xir, yir, interp='nn')
extrapolate_nans(xir,yir,zir)

m, n = 10, 10     # The shape of the window array
# Calculating a couple of smoothed data.
win = numpy.ones((m, n))


####plot
#zis = ndimage.gaussian_filter(zir, sigma=25.0, order=0)
zis = moving_average_2d(zir, win)
xi2= xir.ravel()
yi2= yir.ravel()
zi2= zis.ravel()
zi1= zir.ravel()


fig=plt.figure()
plt.tricontourf(xi2, yi2,zi1,200)
plt.axis((min(xi2),max(xi2),min(yi2),max(yi2)))
plt.title('Before Smothing')
plt.xlabel('X-COORD')
plt.ylabel('Y-COORD')
plt.colorbar()

fig=plt.figure()
plt.tricontourf(xi2, yi2,zi2,200)
plt.axis((min(xi2),max(xi2),min(yi2),max(yi2)))
plt.title('After Smothing')
plt.xlabel('X-COORD')
plt.ylabel('Y-COORD')
plt.colorbar()


fig=plt.figure()
plt.tricontourf(xi2, yi2,zi2,200)
plt.axis((min(xi2),max(xi2),min(yi2),max(yi2)))
plt.plot(x1,y1, 'k. ')
plt.title('After Smothing with original data location')
plt.xlabel('X-COORD')
plt.ylabel('Y-COORD')
plt.colorbar()

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(xi2, yi2, zi2, cmap=cm.jet, linewidth=0.0)
plt.show()