離散的な観測点データから2次元grid作成 綺麗に内挿したい

Question

離散的な観測点データから2次元gridデータへ変換をしたのですが、
もっと見た目を綺麗にしたいです。

試してみたこと
griddata.py で定義された方法で、x,y,z をgrid化しました。
grid の値は、griddata.pyのsizeで定義された範囲で取得される点の中央値（median）を返しています。

#!/usr/bin/env python


import pandas as pd
import griddata
import numpy as np
import sys
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, maskoceans
df = pd.read_csv(sys.argv[1])

#filter
pm10 = df[(df.pm10 != "-")]
# データの準備
X = pm10.lon.values.astype(np.float32)
Y = pm10.lat.values.astype(np.float32)
Z = pm10.pm10.values.astype(np.float32)

print(X.min(),X.max())
print(Y.min(),Y.max())


binsize = 0.01
extent = (124.50,129.50,33.10,39.10)
grid, bins, binloc,xi,yi = griddata.griddata(X, Y, Z, extent,binsize=binsize)  # see this routine's docstring
xx,yy = np.meshgrid(xi,yi)
print(grid.shape)

# minimum values for colorbar. filter our nans which are in the grid
zmin    = grid[np.where(np.isnan(grid) == False)].min()
zmax    = grid[np.where(np.isnan(grid) == False)].max()
grid[np.isnan(grid)] = zmin
grid = maskoceans(xx,yy,grid)
bmap = Basemap(projection='cyl',llcrnrlon=124.50,urcrnrlon=129.50,llcrnrlat=33.10,urcrnrlat=39.10,resolution='i')

bmap.imshow(grid, extent=extent,cmap=plt.cm.rainbow, origin='lower', vmin=0, vmax=zmax, aspect='auto')
bmap.drawcoastlines(linewidth=0.3, linestyle='solid', color='black')
bmap.drawcountries(linewidth=0.3, linestyle='solid', color='black')
bmap.drawparallels(np.arange(-100, 90.0, 1.0), linewidth=0.1, color='black')
bmap.drawmeridians(np.arange(0.0, 360.0, 1.0), linewidth=0.1, color='black')
plt.colorbar() # draw colorbar
plt.scatter(X,Y,marker='o',c='b',s=5)
plt.title('Obs.Analysis TEST grid=0.01')
#plt.show()
plt.savefig("obs_grid_0p01.png",bbox_inches='tight')

# griddata.py - 2010-07-11 ccampo
import numpy as np

def griddata(x, y, z, extent,binsize=0.01, retbin=True, retloc=True):
    """
    Place unevenly spaced 2D data on a grid by 2D binning (nearest
    neighbor interpolation).
    
    Parameters
    ----------
    x : ndarray (1D)
        The idependent data x-axis of the grid.
    y : ndarray (1D)
        The idependent data y-axis of the grid.
    z : ndarray (1D)
        The dependent data in the form z = f(x,y).
    binsize : scalar, optional
        The full width and height of each bin on the grid.  If each
        bin is a cube, then this is the x and y dimension.  This is
        the step in both directions, x and y. Defaults to 0.01.
    retbin : boolean, optional
        Function returns `bins` variable (see below for description)
        if set to True.  Defaults to True.
    retloc : boolean, optional
        Function returns `wherebins` variable (see below for description)
        if set to True.  Defaults to True.
   
    Returns
    -------
    grid : ndarray (2D)
        The evenly gridded data.  The value of each cell is the median
        value of the contents of the bin.
    bins : ndarray (2D)
        A grid the same shape as `grid`, except the value of each cell
        is the number of points in that bin.  Returns only if
        `retbin` is set to True.
    wherebin : list (2D)
        A 2D list the same shape as `grid` and `bins` where each cell
        contains the indicies of `z` which contain the values stored
        in the particular bin.

    Revisions
    ---------
    2010-07-11  ccampo  Initial version
    """
    # get extrema values.
    #xmin, xmax = x.min(), x.max()
    #ymin, ymax = y.min(), y.max()
    xmin,xmax = extent[0],extent[1]
    ymin,ymax = extent[2],extent[3]
    # make coordinate arrays.
    xo      = np.arange(xmin, xmax+binsize, binsize)
    yo      = np.arange(ymin, ymax+binsize, binsize)
    xi, yi  = np.meshgrid(xo,yo)

    # make the grid.
    grid       = np.zeros(xi.shape, dtype=x.dtype)
    nrow, ncol = grid.shape
    if retbin: bins = np.copy(grid)

    # create list in same shape as grid to store indices
    if retloc:
        wherebin = np.copy(grid)
        wherebin = wherebin.tolist()

    # fill in the grid.
    for row in range(nrow):
        for col in range(ncol):
            xc = xi[row, col]    # x coordinate.
            yc = yi[row, col]    # y coordinate.
            
            # find the position that xc and yc correspond to.
            posx = np.abs(x - xc)
            posy = np.abs(y - yc)
            size = 0.5
            #ibin = np.logical_and(posx < binsize/2., posy < binsize/2.)
            ibin = np.logical_and(posx < size/2., posy < size/2.)
            ind  = np.where(ibin == True)[0]
            
            # fill the bin.
            bin = z[ibin]
            if retloc: wherebin[row][col] = ind
            if retbin: bins[row, col] = bin.size
            if bin.size != 0:
                binval         = np.median(bin)
                grid[row, col] = binval
            else:
                grid[row, col] = np.nan   # fill empty bins with nans.

    # return the grid
    if retbin:
        if retloc:
#            return grid, bins, wherebin
            return grid, bins, wherebin,xo,yo
        else:
            return grid, bins
    else:
        if retloc:
            return grid, wherebin
        else:
            return grid

griddata.py size=0.5 の場合

griddata.py size=1.0 の場合

です。
なんとなく良さげにできたんですけど、sizeを大きくしてその中央値をとると
当たり前ですが、平均化されてしまいます。
これに、中央値ではなく、距離による重み付けなどをしてメリハリつけてもっと綺麗に
内挿できればいいと考えています。

どなたかアドバイスいただければ幸いです。

Comments

[hayataka2049]

多重投稿ミス

[hayataka2049]

多重投稿ミス

[hayataka2049]

なにか学術的な用途に使う図なのかなと思いますが、勝手に内挿しちゃっていいんですか？

[shiroikuma]

夏季休暇で返答遅れました。本件は学術的に使うものではなく、面的なデータが必要で、いくつかパターンを作成するのが目的です。内挿方法を数種試してみたい、というのが本件の趣旨です。

Answer 1

使用経験がないので不確かですが、scipy.interpolate.griddataを使うといいかもしれません。

sampleコードを見つけましたので、確認してください

https://scipython.com/book/chapter-8-scipy/examples/two-dimensional-interpolation-with-scipyinterpolategriddata/

Answer 2

https://stackoverflow.com/questions/3104781/inverse-distance-weighted-idw-interpolation-with-python

こちらに投稿するのが遅くなりましたがこちら解決済みです。
調べたところ、IDW法を使い離散的な点を内装する手法がいいとのことで
こちらを採用しました。
他に、rbf,krikingの方法もありましたが、検証した結果、観測値点場所、観測地点数
よりIDW法を採用しました。
ありがとうございました。