Friday, January 15, 2016

Extract Surface at input XY 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
import scipy.ndimage as ndimage

#-- 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('input.txt')
x1 = data[:,0]
y1 = data[:,1]
z1 = data[:,4]

### gridding
numcols, numrows = 100, 100
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 = 15, 15    # 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)

###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)))
xp=numcols*(x1-xmin)/(xmax-xmin)
yp=numrows*(y1-ymin)/(ymax-ymin)

###Extract the values along the line, using cubic interpolation (pixel)
z1s = scipy.ndimage.map_coordinates(zis, np.vstack((xp,yp)))

###save result to file
gabung=np.column_stack((x1,y1,z1,z1s))
np.savetxt('moveave.txt', gabung, fmt='%1.0f')

No comments: