想要从蒙版数组平滑轮廓 [英] want to smooth a contour from a masked array
本文介绍了想要从蒙版数组平滑轮廓的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我有一个被遮罩的数组，matplotlib.plt.contourf使用该数组将温度轮廓投影到glabal地图上。我正在尝试平滑轮廓，但是不幸的是，所有提议的解决方案似乎都无法处理遮罩数组。我测试了以下解决方案：
-。基本上，他使用Inpaint算法对缺失值进行插值，并生成有效的数组进行过滤。
代码在帖子的结尾，您可以将其保存到站点程序包中（或其他方式）并加载为模块（即 inpaint.py ）：
导入油漆
填充=油漆.replace_nans（NANMask，5，0.5，2，'idw'）
我对结果，我想它将适合丢失的温度值。这里还有下一个版本： github 但代码会
为便于参考，易于使用和保存，我将代码（初始版本）放在此处：
/ p>
＃-*-编码：utf-8-*-
各种模块实用程序和帮助器功能
numpy输入为np
#cimport numpy输入为np
#cimport cython
DTYPEf = np.float64
#ctypedef np.float64_t DTYPEf_t
DTYPEi = np.int32
#ctypedef np.int32_t DTYPEi_t
#@cython.boundscheck（False）＃打开整个函数的边界检查
#@cython.wraparound（False）＃整个函数的边界检查
def replace_nans（array，max_iter，tol，kernel_size = 1，method ='localmean'） ：
使用iter替换数组中的NaN元素图像修复算法。
该算法如下：
1）对于输入数组中的每个元素，将其替换为不是NaN本身的相邻元素的加权平均
。权重取决于方法类型的
。如果``method = localmean''权重等于1 /（（（2 * kernel_size + 1）** 2 -1）
2）如果存在相邻的NaN元素，则需要多次迭代。
如果是这种情况，信息会从丢失的
区域的边缘反复传播，直到变化低于某个阈值为止。
参数
----------
数组：2d np.ndarray
包含必须替换的NaN元素的数组
max_iter：int
迭代次数
kernel_size：int
内核大小，默认值为1
方法：str
用于替换无效值的方法。有效选项为
`localmean`，'idw'。
返回
-------
fill：2d np.ndarray
输入数组的副本，其中NaN元素已被替换。
＃cdef int i，j，I，J，it，n，k，l
＃cdef int n_invalids
已填充= np.empty（[array.shape [0]，array.shape [1]]，dtype = DTYPEf）
内核= np.empty（（2 * kernel_size + 1，2 * kernel_size + 1），dtype = DTYPEf）
＃cdef np.ndarray [np.int_t，ndim = 1] inans
＃cdef np.ndarray [np.int_t，ndim = 1] jnans
＃数组为NaN的索引
inans，jnans = np.nonzero（np.isnan（array））
＃NaN元素的数量
n_nans = len（inans ）
＃包含替换值以检查收敛性的数组
replace_new = np.zeros（n_nans，dtype = DTYPEf）
replace_old = np.zeros（n_nans，dtype = DTYPEf ）
＃取决于内核类型，如果方法=='localmean'，则填充内核数组
：
print'kernel_size'，kernel_size
for i范围（2 * kernel_size + 1）中的$：
for范围（2 * kernel_size + 1）中的j：
kernel [i，j] = 1
打印内核，内核
elif方法=='idw'：
内核= np.array（[[0，0.5，0.5，0.5,0]，
[0.5,0.75,0.75,0.75,0.5]，
[0.5,0.75,1,0.75,0.5]，
[0.5,0.75,0.75,0.5,1]，
[0，0.5，0.5，0.5，0.5，0]]）
打印内核，'kernel'
else：
提高ValueError（'方法无效。应该是localmean之一。）
＃用范围为i的输入元素
填充新数组（array.shape [0]）：范围内的j的
（array.shape [1]）：
fill [i，j] = array [i，j]
＃进行多次通过
＃直到达到收敛
为范围（max_iter）：
打印'iteration'，
＃为每个NaN元素
为范围k中的k（n_nans）：
i = inans [k]
j = jnans [k]
＃初始化为零
fill [i，j] = 0.0
n = 0
＃在I在范围（2 * kernel_size + 1）中的内核
：J在范围（2 * kernel_size + 1）中的
：
＃如果我们不在边界之外如果i + I-kernel_size<
array.shape [0]和i + I-kernel_size> = 0：如果j + J-kernel_size< array.shape [1]和j + J-kernel_size> = 0：
＃如果邻居元素本身不是NaN。
如果已填充[i + I内核大小，j + J内核大小] ==已填充[i + I内核大小，j + J内核大小]：
＃不求和
如果I-kernel_size！= 0和J-kernel_size！= 0：
＃卷积具有原始数组
的内核fill [i，j] = filled [i，j] + fill [i + I-kernel_size，j + J-kernel_size] * kernel [I，J]
n = n + 1 * kernel [I，J]
＃将值除以有效数n = 0时添加元素
的数量：
fill [i，j] =填充[i，j] / n
replace_new [k] =填充[i，j]
其他：
fill [i，j] = np.nan
＃检查替换后的
＃的元素之间的均方差是否在某个公差
以下打印'tolerance'，np.mean（（replaced_new-re place_old）** 2）
如果np.mean（（replaced_new-replaced_old）** 2）< tol：
中断
其他：
对于范围内的l（n_nans）：
replace_old [l] = replace_new [l]
返回已填充
def sincinterp（image，x，y，kernel_size = 3）：
在像素之间的中间位置重新采样图像。
此函数使用基本插值公式，它限制了
在重新采样过程中的信息丢失，它使用了数量有限的
个相邻像素。
新的图像：math：ʻim ^ +`在小数位置： math：`x`和：math：`y`的计算公式为：
.. math ::
im ^ +（x，y）= \sum_ {i = -\mathtt {kernel \_size}} ^ {i = \mathtt {kernel\_size}} \sum_ {j = -\mathtt {kernel\_size}} ^ {j = \mathtt {kernel\_size}} \mathtt {image}（i，j）sin [\pi（i-\mathtt {x}）] sin [\pi（j-\mathtt {y}）] / \pi（i- \mathtt {x}）/ \pi（j-\mathtt {y}）
参数
----------
image：np.ndarray， dtype np.int32
图像数组
x：二维np.ndarray of floats
包含小数像素行
位置进行图像插值的数组
y：二维np.ndarray of floats
包含小数像素列
位置进行图像插值的数组
kernel_size：int
插值是在每个插值点附近的（2 * kernel_size + 1）*（2 * kernel_size + 1）
子矩阵上执行的。
返回
-------
im：np.ndarray，dtype np.float64
image在指定点$ b $的插值b由``x''和``y''
＃索引
＃cdef int i，j，I，J
＃输出数组
r = np.zeros（[x.shape [0]，x.shape [1]]，dtype = DTYPEf）
＃快速pi
pi = 3.1419
＃对于范围（I.x.shape [0]）中的I，输出阵列的每个点
：对于范围（x.shape [1]中的J，
]）：
＃循环遍历所有相邻网格点
for i in range（int（x [I，J]）-kernel_size，int（x [I，J]）+ kernel_size +1）：
对于范围内的j（int（y [I，J]）-kernel_size，int（y [I，J]）+ kernel_size + 1）：
＃检查我们是否在如果i> = 0且i< = image.shape [0]和j> = 0且j< = image.shape [1]，则边界
：
if（ix [I ，J]] == 0.0和（jy [I，J]）== 0.0：
r [I，J] = r [I，J] +图像[i，j]
elif（ix [I，J]）== 0.0：
r [I，J] = r [I ，J] + image [i，j] * np.sin（pi *（jy [I，J]））/（pi *（jy [I，J]））
elif（jy [I，J ]）== 0.0：
r [I，J] = r [I，J] + image [i，j] * np.sin（pi *（ix [I，J]））/（pi *（ ix [I，J]））
else：
r [I，J] = r [I，J] + image [i，j] * np.sin（pi *（ix [I，J ]））* np.sin（pi *（jy [I，J]））/（pi * pi *（ix [I，J]）*（jy [I，J]））
return r
#cdef来自 math.h的外部变量：
＃double sin（double）
I have a masked array which is used by matplotlib.plt.contourf to project a temperature contour on a glabal map. I was trying to smooth the contour, but unfortunately none of the proposed solutions seems to be able to handle masked array. I tested these solutions:
-scipy.ndimage.gaussian_filter - moving averages
- scipy.ndimage.zoom
none of them works(they count in the masked values also). Is there any way I can smooth my contour on maskedArray
I have added this part after trying the proposed 'inpaint' solution and the results were unchanged. here is the code (if it helps)
import Scientific.IO.NetCDF as S
import mpl_toolkits.basemap as bm
import numpy.ma as MA
import numpy as np
import matplotlib.pyplot as plt
import inpaint
def main():
fileobj = S.NetCDFFile('Bias.ANN.tas_A1_1.nc', mode='r')
# take the values
set1 = {'time', 'lat', 'lon'}
set2 = set(fileobj.variables.keys())
set3 = set2 - set1
datadim = set3.pop()
print "******************datadim: "+datadim
data = fileobj.variables[datadim].getValue()[0,:,:]
lon = fileobj.variables['lon'].getValue()
lat = fileobj.variables['lat'].getValue()
fileobj.close()
data, lon = bm.shiftgrid(180.,data, lon,start=False)
data = MA.masked_equal(data, 1.0e20)
#data2 = inpaint.replace_nans(data, 10, 0.25, 2, 'idw')
#- Make 2-D longitude and latitude arrays:
[lon2d, lat2d] =np.meshgrid(lon, lat)
#- Set up map:
mapproj = bm.Basemap(projection='cyl',
llcrnrlat=-90.0, llcrnrlon=-180.00,
urcrnrlat=90.0, urcrnrlon=180.0)
mapproj.drawcoastlines(linewidth=.5)
mapproj.drawmapboundary(fill_color='.8')
#mapproj.drawparallels(N.array([-90, -45, 0, 45, 90]), labels=[1,0,0,0])
#mapproj.drawmeridians(N.array([0, 90, 180, 270, 360]), labels=[0,0,0,1])
lonall, latall = mapproj(lon2d, lat2d)
cmap=plt.cm.Spectral
#- Make a contour plot of the temperature:
mymapf = plt.contourf(lonall, latall, data, 20, cmap=cmap)
#plt.clabel(mymapf, fontsize=12)
plt.title(cmap.name)
plt.colorbar(mymapf, orientation='horizontal')
plt.savefig('sample2.png', dpi=150, edgecolor='red', format='png', bbox_inches='tight', pad_inches=.2)
plt.close()
if __name__ == "__main__":
main()
I am comparing the output from this code (the first figure), with output of the same datafile from Panoply. Zoomin in and looking more precisely it seems like it is not the smoothness problem, but the pyplot model provides one stripe slimmer, or the contours are cut earlier (the outer boundaries shows this clearly, and inner contours are different due to this fact). It makes it to look like that the pyplot model is not as smooth as the Panoply one. how can I get (nearly) the same model? Am I distinguishing it right?
解决方案
I had similar problem and google pointed me to this: blog post. Basically he's using inpaint algorithm to interpolate missing values and produce valid array for filtering.
The code is at the end of the post, and you can save it to site-packages (or else) and load it as module (i.e. inpaint.py):
import inpaint
filled = inpaint.replace_nans(NANMask, 5, 0.5, 2, 'idw')
I'm happy with the result, and I guess it will suite missing temperature values just fine. There is also next version here: github but code will need some cleaning for general usage as it's part of a project.
For reference, easy use and preservation sake I'll post the code (of initial version) here:
# -*- coding: utf-8 -*-
"""A module for various utilities and helper functions"""
import numpy as np
#cimport numpy as np
#cimport cython
DTYPEf = np.float64
#ctypedef np.float64_t DTYPEf_t
DTYPEi = np.int32
#ctypedef np.int32_t DTYPEi_t
#@cython.boundscheck(False) # turn of bounds-checking for entire function
#@cython.wraparound(False) # turn of bounds-checking for entire function
def replace_nans(array, max_iter, tol,kernel_size=1,method='localmean'):
"""Replace NaN elements in an array using an iterative image inpainting algorithm.
The algorithm is the following:
1) For each element in the input array, replace it by a weighted average
of the neighbouring elements which are not NaN themselves. The weights depends
of the method type. If ``method=localmean`` weight are equal to 1/( (2*kernel_size+1)**2 -1 )
2) Several iterations are needed if there are adjacent NaN elements.
If this is the case, information is "spread" from the edges of the missing
regions iteratively, until the variation is below a certain threshold.
Parameters
----------
array : 2d np.ndarray
an array containing NaN elements that have to be replaced
max_iter : int
the number of iterations
kernel_size : int
the size of the kernel, default is 1
method : str
the method used to replace invalid values. Valid options are
`localmean`, 'idw'.
Returns
-------
filled : 2d np.ndarray
a copy of the input array, where NaN elements have been replaced.
"""
# cdef int i, j, I, J, it, n, k, l
# cdef int n_invalids
filled = np.empty( [array.shape[0], array.shape[1]], dtype=DTYPEf)
kernel = np.empty( (2*kernel_size+1, 2*kernel_size+1), dtype=DTYPEf )
# cdef np.ndarray[np.int_t, ndim=1] inans
# cdef np.ndarray[np.int_t, ndim=1] jnans
# indices where array is NaN
inans, jnans = np.nonzero( np.isnan(array) )
# number of NaN elements
n_nans = len(inans)
# arrays which contain replaced values to check for convergence
replaced_new = np.zeros( n_nans, dtype=DTYPEf)
replaced_old = np.zeros( n_nans, dtype=DTYPEf)
# depending on kernel type, fill kernel array
if method == 'localmean':
print 'kernel_size', kernel_size
for i in range(2*kernel_size+1):
for j in range(2*kernel_size+1):
kernel[i,j] = 1
print kernel, 'kernel'
elif method == 'idw':
kernel = np.array([[0, 0.5, 0.5, 0.5,0],
[0.5,0.75,0.75,0.75,0.5],
[0.5,0.75,1,0.75,0.5],
[0.5,0.75,0.75,0.5,1],
[0, 0.5, 0.5 ,0.5 ,0]])
print kernel, 'kernel'
else:
raise ValueError( 'method not valid. Should be one of `localmean`.')
# fill new array with input elements
for i in range(array.shape[0]):
for j in range(array.shape[1]):
filled[i,j] = array[i,j]
# make several passes
# until we reach convergence
for it in range(max_iter):
print 'iteration', it
# for each NaN element
for k in range(n_nans):
i = inans[k]
j = jnans[k]
# initialize to zero
filled[i,j] = 0.0
n = 0
# loop over the kernel
for I in range(2*kernel_size+1):
for J in range(2*kernel_size+1):
# if we are not out of the boundaries
if i+I-kernel_size < array.shape[0] and i+I-kernel_size >= 0:
if j+J-kernel_size < array.shape[1] and j+J-kernel_size >= 0:
# if the neighbour element is not NaN itself.
if filled[i+I-kernel_size, j+J-kernel_size] == filled[i+I-kernel_size, j+J-kernel_size] :
# do not sum itself
if I-kernel_size != 0 and J-kernel_size != 0:
# convolve kernel with original array
filled[i,j] = filled[i,j] + filled[i+I-kernel_size, j+J-kernel_size]*kernel[I, J]
n = n + 1*kernel[I,J]
# divide value by effective number of added elements
if n != 0:
filled[i,j] = filled[i,j] / n
replaced_new[k] = filled[i,j]
else:
filled[i,j] = np.nan
# check if mean square difference between values of replaced
#elements is below a certain tolerance
print 'tolerance', np.mean( (replaced_new-replaced_old)**2 )
if np.mean( (replaced_new-replaced_old)**2 ) < tol:
break
else:
for l in range(n_nans):
replaced_old[l] = replaced_new[l]
return filled
def sincinterp(image, x, y, kernel_size=3 ):
"""Re-sample an image at intermediate positions between pixels.
This function uses a cardinal interpolation formula which limits
the loss of information in the resampling process. It uses a limited
number of neighbouring pixels.
The new image :math:`im^+` at fractional locations :math:`x` and :math:`y` is computed as:
.. math::
im^+(x,y) = \sum_{i=-\mathtt{kernel\_size}}^{i=\mathtt{kernel\_size}} \sum_{j=-\mathtt{kernel\_size}}^{j=\mathtt{kernel\_size}} \mathtt{image}(i,j) sin[\pi(i-\mathtt{x})] sin[\pi(j-\mathtt{y})] / \pi(i-\mathtt{x}) / \pi(j-\mathtt{y})
Parameters
----------
image : np.ndarray, dtype np.int32
the image array.
x : two dimensions np.ndarray of floats
an array containing fractional pixel row
positions at which to interpolate the image
y : two dimensions np.ndarray of floats
an array containing fractional pixel column
positions at which to interpolate the image
kernel_size : int
interpolation is performed over a ``(2*kernel_size+1)*(2*kernel_size+1)``
submatrix in the neighbourhood of each interpolation point.
Returns
-------
im : np.ndarray, dtype np.float64
the interpolated value of ``image`` at the points specified
by ``x`` and ``y``
"""
# indices
# cdef int i, j, I, J
# the output array
r = np.zeros( [x.shape[0], x.shape[1]], dtype=DTYPEf)
# fast pi
pi = 3.1419
# for each point of the output array
for I in range(x.shape[0]):
for J in range(x.shape[1]):
#loop over all neighbouring grid points
for i in range( int(x[I,J])-kernel_size, int(x[I,J])+kernel_size+1 ):
for j in range( int(y[I,J])-kernel_size, int(y[I,J])+kernel_size+1 ):
# check that we are in the boundaries
if i >= 0 and i <= image.shape[0] and j >= 0 and j <= image.shape[1]:
if (i-x[I,J]) == 0.0 and (j-y[I,J]) == 0.0:
r[I,J] = r[I,J] + image[i,j]
elif (i-x[I,J]) == 0.0:
r[I,J] = r[I,J] + image[i,j] * np.sin( pi*(j-y[I,J]) )/( pi*(j-y[I,J]) )
elif (j-y[I,J]) == 0.0:
r[I,J] = r[I,J] + image[i,j] * np.sin( pi*(i-x[I,J]) )/( pi*(i-x[I,J]) )
else:
r[I,J] = r[I,J] + image[i,j] * np.sin( pi*(i-x[I,J]) )*np.sin( pi*(j-y[I,J]) )/( pi*pi*(i-x[I,J])*(j-y[I,J]))
return r
#cdef extern from "math.h":
# double sin(double)
这篇关于想要从蒙版数组平滑轮廓的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
