Friday, January 15, 2016

2D Smoothing with Moving Average and Save Result to file 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

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), 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)

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

####plot
#zis = ndimage.gaussian_filter(zi, sigma=25.0, order=0)
zis = moving_average_2d(zi, win)
#print xi.shape
#print yi.shape
#print zis.shape

###save result to file
gabung=np.column_stack((xi.reshape((numel,1)),yi.reshape((numel,1)),zis.reshape((numel,1))))
np.savetxt('moveave.txt', gabung, fmt='%1.0f')
fig=plt.figure()
ax=fig.add_subplot(1,2,1)
plt.contour(xi, yi, zi, 20, linewidths = 0.5, colors = 'k')
plt.contourf(xi,yi,zis)
ax=fig.add_subplot(1,2,2)
plt.contour(xi, yi, zis, 20, linewidths = 0.5, colors = 'k')
plt.contourf(xi,yi,zis)
plt.colorbar()
plt.show()

No comments: