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()
No comments:
Post a Comment