Python basemap.interp函数代码示例

    def plot(self):

        self.data, lonwrap = addcyclic(self.data, self.lons)

        # Sort latitudes and data
        lat_idx = np.argsort(self.lats)
        self.lats = self.lats[lat_idx]
        self.data = self.data[lat_idx]

        data_lon_min = min(lonwrap)
        data_lon_max = max(lonwrap)
        data_lat_min = min(self.lats)
        data_lat_max = max(self.lats)

        new_lons = np.arange(data_lon_min - 1.0, data_lon_max + 1.0, 1.0)
        new_lats = np.arange(data_lat_min - 1.0, data_lat_max + 1.0, 1.0)

        x, y = self.m(*np.meshgrid(new_lons[:], new_lats[:]))

        # Two pass interpolation to deal with the mask.
        # First pass does bilinear, the next does nearest neighbour
        # interpolation.
        # It's not clear this is working, and the problem is likely
        # solved by ensuring the right mask is used!
        data_bl = interp(self.data, lonwrap[:], self.lats[:], x, y, order=1)
        data_nn = interp(self.data, lonwrap[:], self.lats[:], x, y, order=0)

        data_bl[data_nn.mask == 1] = data_nn[data_nn.mask == 1]

        if self.parameters.has_key("color_levels"):
            self.__print_custom_color_plot(x, y, data_bl)
        else:
            self.__print_cmap_plot(x, y, data_bl)

        return self.main_render

    def extract_transect(self, dataout_line, data=None, data_x=None, data_y=None):
        if data is 'depth':
            data=self.depth
        elif data is 'strike':
            data=self.strike
        elif data is 'dip':
            data=self.dip
        elif data is None:
            print("Error in extract transect: need to specify input for 'data'")
        else:
            data=data

        if data_x is not None:
            data_x=data_x
        else:
            data_x=self.x
        if data_y is not None:
            data_y=data_y
        else:
            data_y=self.y

        dataout1 = interp(data, data_x, data_y, dataout_line[:,0],dataout_line[:,1], order=1)
        dataout2 = interp(data, data_x, data_y, dataout_line[:,0],dataout_line[:,1], order=0)

        for i in range(0,np.size(dataout1)):
            if dataout1[i] is np.ma.masked:
                if dataout2[i] is not np.ma.masked:
                    dataout1[i] = dataout2[i]
                else:
                    r = i

                    while dataout2[r] is np.ma.masked:
                        if r < np.size(dataout1) - 1:
                            r += 1
                    try:
                        right = dataout2[r]
                    except IndexError:
                        pass



                    l = i
                    while dataout2[l-1] is np.ma.masked:
                        l += -1
                    try:
                        left = dataout2[l-1]
                    except IndexError:
                        pass

                    dataout1[i] = np.average([right,left])

        return dataout1

def regrid(x, y, arr, inc_by=2):
    """Regrid a 2d array increasing its resolution."""
    ny, nx = arr.shape
    xi = np.linspace(x.min(), x.max(), inc_by * len(x))
    yi = np.linspace(y.min(), y.max(), inc_by * len(y))
    xx, yy = np.meshgrid(xi, yi)
    arr = np.ma.masked_invalid(arr)
    arr1 = bm.interp(arr, x, y, xx, yy, order=0) # nearest neighb.
    arr2 = bm.interp(arr, x, y, xx, yy, order=1) # linear interp.
    ind = np.where(arr2 == 0)                    #<<<<< check!
    try:
        arr2[ind] = arr1[ind]
    except:
        pass
    return [xi, yi, arr2]

def interpdata(data, lats, lons):

    '''
    Interpola dados para 1 grau
    Os dados de entrada devem ter 3 dimenões: tempo, lat, lon

    :param: data - Dados com 3 dimensões
    :type param: numpy array
    :param: lats - Latitudes a serem interpoladas
    :type param: numpy array 1d
    :param: lons - Longitudes a serem interpoladas
    :type param: numpy array 1d
    '''

    # Criando grade de 1 grau
    newlats = np.linspace(-90, 90, 181)
    newlons = np.linspace(-180, 179, 360)
    x, y = np.meshgrid(newlons, newlats)

    # Interpola dados
    newdata = np.empty((int(data.shape[0]), int(len(newlats)), int(len(newlons))))
    for i in range(0, int(data.shape[0])):
        newdata[i, :, :] = interp(data[i, :, :], lons, lats, x, y, order=1)

    return newdata, newlats, newlons

def interpolate(data,navlon,navlat,interp=None):
    """Perform a spatial interpolation if required; return x_reg,y_reg,data_reg.
    """
    e1,e2 = e1e2(navlon,navlat) # ideally we would like e1u and not e1t...
    x1d_in = e1[0,:].cumsum() - e1[0,0]
    y1d_in = e2[:,0].cumsum() - e2[0,0]
    x2d_in,y2d_in = npy.meshgrid(x1d_in,y1d_in)
    # print x1d_in
    if interp is None or interp=='0':
       return x2d_in, y2d_in, data.copy()
    elif interp=='basemap': # only for rectangular grid...
       from mpl_toolkits import basemap
       x1d_reg=npy.linspace(x1d_in[0],x1d_in[-1],len(x1d_in))
       y1d_reg=npy.linspace(y1d_in[0],y1d_in[-1],len(y1d_in))
       x2d_reg,y2d_reg = npy.meshgrid(x1d_reg,y1d_reg)
       data_reg=basemap.interp(data,x1d_in,y1d_in,x2d_reg,y2d_reg,checkbounds=False,order=1)
       return x2d_reg,y2d_reg,data_reg
    elif interp=='scipy': # only for rectangular grid...
       import scipy.interpolate
       x1d_reg=npy.linspace(x1d_in[0],x1d_in[-1],len(x1d_in))
       y1d_reg=npy.linspace(y1d_in[0],y1d_in[-1],len(y1d_in))
       x2d_reg,y2d_reg = npy.meshgrid(x1d_reg,y1d_reg)
       interp = scipy.interpolate.interp2d(x1d_in, y1d_in,data, kind='linear')
       a1d = interp(x2d_reg[0,:],y2d_reg[:,0])
       data_reg = npy.reshape(a1d,y2d_reg.shape)
       #test_plot(x2d_in,y2d_in,data)
       #test_plot(x2d_reg,y2d_reg,data_reg)
       return x2d_reg,y2d_reg,data_reg

def regrid(lon, lat, dhdt, factor=10):
    m, n = len(lon), len(lat)
    lon2 = np.linspace(lon[0], lon[-1], m*factor)
    lat2 = np.linspace(lat[0], lat[-1], n*factor)
    xx, yy = np.meshgrid(lon2, lat2)
    dhdt2 = interp(dhdt, lon, lat, xx, yy, order=1)                          # good!!!
    return lon2, lat2, dhdt2

def extened_grid(zi,x1,y1,zoom=2):
    '''
    xinterval : X插值的间隔
    yinterval : Y 插值的间隔
    扩展网格区域zoom为扩展倍数
    '''
    #print(x1)
    nx = np.size(x1)
    ny = np.size(y1)
    x2 = np.linspace(x1.min(), x1.max(), nx * zoom)
    y2 = np.linspace(y1.min(), y1.max(), ny * zoom)
    xi,yi = np.meshgrid(x2,y2)

    #插值方法1 Zoom方法
    #from scipy import ndimage
    #z2 = ndimage.interpolation.zoom(zi[:,:], zoom)

    #插值方法2 basemap.interp方法
    from mpl_toolkits.basemap import interp
    z2 = interp(zi, x1, y1, xi, yi, checkbounds=True, masked=False, order=1)

    #插值方法3 interpolate.RectBivariateSpline 矩形网格上的样条逼近。
    # Bivariate spline approximation over a rectangular mesh
    #from scipy import interpolate
    #sp = interpolate.RectBivariateSpline(y1,x1,zi,kx=1, ky=1, s=0)
    #z2 = sp(y2,x2)

    #sp = interpolate.LSQBivariateSpline(y1,x1,zi)
    #z2 = sp(y2,x2)

    #terpolate.LSQBivariateSpline?

    print('extend shapes:=',z2.shape,xi.shape,yi.shape)
    return z2,xi,yi,x2,y2,nx*zoom,ny*zoom

 def extend_interp(datafield):
   # add masked values at southernmost end
   southernlimitmask = ma.masked_all(len(self.olon))
   olat_ext          = np.append(-82.1,self.olat)
   dfield_ext = ma.concatenate([ma.column_stack(southernlimitmask), datafield], 0)
   # f = interp2d(self.olon, olat_ext, dfield_ext)
   # return f(self.pismlon, self.pismlat)
   return interp(dfield_ext, self.olon, olat_ext, self.pismlon, self.pismlat)

 def __call__(self,array):
     masked = ma.is_masked(array)
     if self.method is 'basemap':
        return basemap.interp(array, self.xin, self.yin, self.xout, self.yout, checkbounds=False, masked=masked, order=1)
     elif self.method is 'scipy':
        import scipy.interpolate
        interp = scipy.interpolate.interp2d(self.xin, self.yin, array, kind='linear')
        a1d = interp(self.xout[0,:],self.yout[:,0])
        return npy.reshape(a1d,self.yout.shape)

def resample_slice(slice_, grid_lon, grid_lat, order=1):
    """
    Resample a single time slice of a larger xr.DataArray

    :param slice: xr.DataArray single slice
    :param grid_lon: meshgrid of longitudes for the new grid
    :param grid_lat: meshgrid of latitudes for the new grid
    :param order: Interpolation method 0 - nearest neighbour, 1 - bilinear (default), 3 - cubic spline
    :return: xr.DataArray, resampled slice
    """
    result = basemap.interp(slice_.values, slice_['lon'].data, slice_['lat'].data, grid_lon, grid_lat)
    return xr.DataArray(result)

def interp_CRU(path, fname, long_new, lat_new, zip=True, dtype=None):
    """
    Extracts from a CRU file, interpolates it to a non-grid point set.
    """
    from mpl_toolkits import basemap
    long_old, lat_old, data = CRU_extract(path, fname, zip, dtype)
    N_new = len(long_new)
    out_vals = zeros(N_new, dtype=float)

    for i in xrange(N_new):
        out_vals[i] = basemap.interp(data,long_old,lat_old,long_new[i],lat_new[i],order=1)
    return out_vals

def regrid2d(arr3d, x, y, inc_by=2):
    """Regrid 2d time series (3d array) increasing resolution."""
    nt, ny, nx = arr3d.shape
    out = np.empty((nt, inc_by * ny, inc_by * nx), 'f8')
    xi = np.linspace(x.min(), x.max(), inc_by * len(x))
    yi = np.linspace(y.min(), y.max(), inc_by * len(y))
    xx, yy = np.meshgrid(xi, yi)
    arr3d = np.ma.masked_invalid(arr3d)
    for k, field in enumerate(arr3d):
        field1 = bm.interp(field, x, y, xx, yy, order=0) # nearest neighb.
        field2 = bm.interp(field, x, y, xx, yy, order=1) # linear
        ##########################################################
        # soemthing "wierd" when the field is zero
        ind = np.where(field2 == 0) #<<<<< check!
        try:
            field2[ind] = field1[ind]
        except:
            pass
        ##########################################################
        out[k] = field2
    return [out, xi, yi]

def griddata_nearest(x, y, z, xi, yi):
    x = x.astype(np.float32)
    y = y.astype(np.float32)
    z = z.astype(np.float32)
    xi = xi.astype(np.float32)
    yi = yi.astype(np.float32)

    (nx,ny)=xi.shape
    xi, yi = xi.flatten(), yi.flatten()
    from scipy.interpolate import griddata
    interp = griddata((x, y), z,(xi,yi), method='nearest')#linear
    zi = np.reshape(interp(xi, yi),(nx,ny))
    zi = zi.astype(np.float32)
    return zi

def griddata_linear_rbf2(x, y, z, xi, yi,function='linear'):

    x = x.astype(np.float32)
    y = y.astype(np.float32)
    z = z.astype(np.float32)
    xi = xi.astype(np.float32)
    yi = yi.astype(np.float32)

    (nx,ny)=xi.shape
    xi, yi = xi.flatten(), yi.flatten()
    from scipy.interpolate import Rbf
    interp = Rbf(x, y, z, epsilon=1)#linear
    zi = np.reshape(interp(xi, yi),(nx,ny))
    zi = zi.astype(np.float32)
    return zi

def interpolate(tecmap):
    """Interpolate TEC Map."""
    lat2 = np.linspace(tecmap.lat[0][0], tecmap.lat[-1][0],
                       tecmap.lat.shape[0] * 10)
    lon2 = np.linspace(tecmap.lon[0][0], tecmap.lon[0][-1],
                       tecmap.lon.shape[1] * 20)
    lon_inter, lat_inter = np.meshgrid(lon2, lat2)
    tecmap_inter = interp(
        tecmap.value,
        tecmap.lon[0],
        np.flipud(tecmap.lat[:, 0]),
        lon_inter,
        np.flipud(lat_inter),
        order=1)
    return lon_inter, lat_inter, tecmap_inter