griddata

mipylib.numeric.interpolate.griddata(points, values, xi=None, **kwargs)

Interpolate scattered data to grid data.

Parameters:
  • points – (list) The list contains x and y coordinate arrays of the scattered data.

  • values – (array_like) The scattered data array.

  • xi – (list) The list contains x and y coordinate arrays of the grid data. Default is None, the grid x and y coordinate size were both 500.

  • method – (string) The interpolation method. [idw | cressman | nearest | inside_mean | inside_min | inside_max | inside_sum | inside_count | surface | barnes | kriging]

  • fill_value – (float) Fill value, Default is nan.

  • pointnum – (int) Only used for ‘idw’ method. The number of the points to be used for each grid value interpolation.

  • radius – (float) Used for ‘idw’, ‘cressman’ and ‘neareast’ methods. The searching raduis. Default is None in ‘idw’ method, means no raduis was used. Default is [10, 7, 4, 2, 1] in cressman method.

  • convexhull – (boolean) If the convexhull will be used to mask result grid data. Default is False.

Returns:

(array) Interpolated grid data (2-D array)

IDW method example

f = addfile('D:/temp/nc/out.20140421_20140421_JL3KMmeic.nc')
data = f['PM25']
data = data[15,1,:,:]
lon = f['lon'][:,:]
lat = f['lat'][:,:]
#Interpolate data to grid
lon1 = linspace(lon.min(), lon.max(), lon.dimlen(1)*5)
lat1 = linspace(lat.min(), lat.max(), lat.dimlen(0)*5)
data1 = griddata((lon, lat), data, xi=(lon1, lat1), method='idw', convexhull=True)[0]
lon_g,lat_g = meshgrid(lon1, lat1)
#Plot
axesm()
mlayer = shaperead('D:/temp/map/jilin.shp')
geoshow(mlayer, edgecolor='r', size=2)
layer = contourfm(lon1, lat1, data1, 20)
scatterm(lon, lat, data, 20, fill=False)
colorbar(layer)
xlim(126.25,126.85)
ylim(43.55,44.15)
grid(True)
title('PM2.5 concentration')
../../../../_images/griddata_convexhull.png

Kriging method example

fn = os.path.join(migl.get_sample_folder(), 'MICAPS', '10101414.000')
f = addfile_micaps(fn)
pr = f['Precipitation6h'][:]
lon = f['Longitude'][:]
lat = f['Latitude'][:]

mask_rect = geolib.polygon([72,72,136,136], [16,55,55,16])
pr,lon,lat = geolib.rmaskout(pr, lon, lat, mask_rect)
idx = where(pr != nan)
pr = pr[idx]
lon = lon[idx]
lat = lat[idx]

#griddata function - interpolate
x = arange(75, 135, 0.5)
y = arange(18, 55, 0.5)
prg = griddata((lon, lat), pr, xi=(x, y), method='kriging')[0]

#Plot
figure(figsize=[700,550], newfig=False)
proj = projinfo(proj='lcc', lon_0=105, lat_1=25, lat_2=47)
axesm(projinfo=proj, position=[0.01, 0.01, 0.99, 0.99], axison=False, gridlabel=False, frameon=False)
geoshow('cn_province', edgecolor='lightgray')
bou1_layer = geoshow('cn_border', facecolor=(0,0,255))
city_layer = geoshow('cn_cities', facecolor='r', size=4)
city_layer.addlabels('NAME', fontname=u'楷体', fontsize=16, yoffset=15, avoidcoll=False)
china_layer = geoshow('china', visible=False)
city_layer.movelabel(u'西宁', -15)
city_layer.movelabel(u'海口', -20, -10)
city_layer.movelabel(u'澳门', 0, -25)
city_layer.movelabel(u'香港', 20, -10)
city_layer.movelabel(u'福州', -10)
city_layer.movelabel(u'合肥', -18)
city_layer.movelabel(u'杭州', 0, -20)
city_layer.movelabel(u'上海', 18)
city_layer.movelabel(u'太原', 0, -20)
city_layer.movelabel(u'天津', 15)
city_layer.movelabel(u'石家庄', -10)
levs = [0.1, 1, 2, 5, 10, 20, 25, 50, 100]
cols = [(255,255,255),(170,240,255),(120,230,240),(200,220,50),(240,220,20),(255,120,10),(255,90,10), \
    (240,40,0),(180,10,0),(120,10,0)]
layer = contourf(x, y, prg, levs, colors=cols)
masklayer(china_layer, [layer])
legend(layer, loc='lower left', frameon=True, title=u'降雨量(mm)', titlefontname=u'黑体')
axism([79, 128, 14, 53])
text(95, 53, u'全国降水量实况图', fontname=u'黑体', fontsize=18)
text(95, 51, u'(2010-10-14 08:00 至 2010-10-14 14:00)', fontname=u'黑体', fontsize=16)

#Add south China Sea
sc_layer = bou1_layer.clone()
axesm(position=[0.15,0.05,0.15,0.2], axison=False, frameon=True)
geoshow(sc_layer, facecolor=(0,0,255))
xlim(106, 123)
ylim(2, 23)
../../../../_images/kriging_griddata.PNG