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

 def _interpolate_structured(self, data, src_lat, src_lon,
                             trg_lat, trg_lon, order=0):
     """
     Interpolate structured data with basemap.interp function.
     """
     reshaped_data = data.reshape((-1, data.shape[-2], data.shape[-1]))
     remapped_data = np.zeros(
         (reshaped_data.shape[0], trg_lat.shape[-2], trg_lat.shape[-1]))
     for i in range(reshaped_data.shape[0]):
         sliced_array = reshaped_data[i, :, :]
         remapped_data[i, :, :] = interp(sliced_array.T, src_lat, src_lon,
                                         trg_lat, trg_lon, order=order)
     remapped_shape = list(data.shape[:-2])+list(remapped_data.shape[-2:])
     remapped_data = remapped_data.reshape(remapped_shape)
     remapped_data = np.atleast_2d(remapped_data)
     return remapped_data

def reproject_data(location, varname, map, lonname='lon', latname='lat', step=1, xsize=100, ysize=100, filter=np.nan):
    nc = Dataset(location)
    latvar = nc.variables[latname]
    lonvar = nc.variables[lonname]
    datavar = nc.variables[varname]

    lons = lonvar[::step]
    lats = latvar[::step]
    if len(datavar.dimensions) == 2:
        data = datavar[::step, ::step]
    elif len(datavar.dimensions) == 3:
        data = datavar[0,::step, ::step]
    
    # Set masked (i.e. land) data to 0.
    # plot_surface ignores masks, and if we set it to NaN, it screws up the colour map
    # TODO: try this again with a custom colour map...
    if filter is not None:
        data[np.where(np.ma.getmask(data) == True)] = filter

    # Now fix the longitude wrapping so that all values go from -180:180
    wrapindex = None
    for i, lon in enumerate(lons):
        if lon > 180:
            lons[i] -= 360
            if wrapindex is None:
                wrapindex = i
    if wrapindex is not None: 
        lons = np.hstack((lons[wrapindex:],lons[:wrapindex]))
        data = np.hstack((data[:,wrapindex:],data[:,:wrapindex]))

    lons_proj, lats_proj = map.makegrid(xsize, ysize)
    
    data_proj = interp(data, lons, lats, lons_proj, lats_proj, checkbounds=False, masked=False, order=1)

    XX, YY = np.meshgrid(np.arange(xsize), np.arange(ysize))
    
    nc.close()
    return (lons_proj, lats_proj, XX, YY, data_proj)

def griddata_scipy_idw(x, y, z, xi, yi,function='linear'):
    '''
    scipy反向距离加权插值
    'multiquadric': sqrt((r/self.epsilon)**2 + 1)  #不能
    'inverse': 1.0/sqrt((r/self.epsilon)**2 + 1) #不能
    'gaussian': exp(-(r/self.epsilon)**2) 不能用来插值
    'linear': r  #能
    'cubic': r**3 #能
    'quintic': r**5  #效果差，勉强能
    'thin_plate': r**2 * log(r)  能可以用用来插值
    '''
    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, function=function,epsilon=2)#linear
    zi = np.reshape(interp(xi, yi),(nx,ny))
    zi = zi.astype(np.float32)
    return zi

def createmap(data, lats, lons, make_edges=False, GC_shift=True,
              vmin=None, vmax=None, latlon=True,
              region=__GLOBALREGION__, aus=False, linear=False,
              clabel=None, colorbar=True, cbarfmt=None, cbarxtickrot=None,
              ticks=None, cbarorient='bottom',
              xticklabels=None,
              set_bad=None, set_under=None, set_over=None,
              pname=None,title=None,suptitle=None, smoothed=False,
              cmapname=None):
    '''
        Pass in data[lat,lon], lats[lat], lons[lon]
        arguments:
            set_bad='blue' #should mask nans as blue
            GC_shift=True #will shift plot half a square left and down
        Returns map, cs, cb
    '''

    # Create a basemap map with region as inputted
    if aus: region=__AUSREGION__
    if __VERBOSE__:
        print("createmap called over %s (S,W,N,E)"%str(region))
        #print("Data %s, %d lats and %d lons"%(str(data.shape),len(lats), len(lons)))

    # First reduce data,lats,lons to the desired region (should save plotting time)
    regionplus=np.array(region) + np.array([-5,-10,5,10]) # add a little padding so edges aren't lost
    lati,loni=util.lat_lon_range(lats,lons,regionplus)
    data=data[lati,:]
    data=data[:,loni]
    lats=lats[lati]
    #print(lons)
    #print(loni)

    lons=lons[loni]


    lllat=region[0]; urlat=region[2]; lllon=region[1]; urlon=region[3]
    m=Basemap(llcrnrlat=lllat, urcrnrlat=urlat, llcrnrlon=lllon, urcrnrlon=urlon,
              resolution='i', projection='merc')

    # plt.colormesh arguments will be added to dictionary
    pcmeshargs={}

    if not linear:
        if __VERBOSE__:
            print('removing %d negative datapoints in createmap'%np.nansum(data<0))
        # ignore warnings of NaN comparison
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore",category =RuntimeWarning)
            data[data<=0] = np.NaN
        pcmeshargs['norm']=LogNorm()

    # Set vmin and vmax if necessary
    if vmin is None:
        vmin=1.05*np.nanmin(data)
    if vmax is None:
        vmax=0.95*np.nanmax(data)

    ## basemap pcolormesh uses data edges
    ##
    lats_e,lons_e=lats,lons
    lats_m,lons_m=lats,lons
    if make_edges:
        if __VERBOSE__: print("Making edges from lat/lon mids")
        nlat,nlon=len(lats), len(lons)
        lats_e=regularbounds(lats)
        lons_e=regularbounds(lons)
        assert nlat == len(lats_e)-1, "regularbounds failed: %d -> %d"%(nlat, len(lats_e))
        assert nlon == len(lons_e)-1, "regularbounds failed: %d -> %d"%(nlon, len(lons_e))
        ## midpoints, derive simply from edges
        lons_m=(lats_e[0:-1] + lats_e[1:])/2.0
        lats_m=(lons_e[0:-1] + lons_e[1:])/2.0
    elif GC_shift: # non edge-based grids need to be shifted left and down by half a box
        latres=lats[3]-lats[2]
        lonres=lons[3]-lons[2]
        lats=lats-latres/2.0
        lons=lons-lonres/2.0
        lats[lats < -89.9] = -89.9
        lats[lats > 89.9]  =  89.9
        lats_e,lons_e=lats,lons
        lats_m,lons_m=lats,lons

    ## interpolate for smoothed output if desired
    ##
    if smoothed:
        factor=5
        if __VERBOSE__: print("Smoothing data, by factor of %d"%factor)
        # 'increase' resolution
        nlats = factor*data.shape[0]
        nlons = factor*data.shape[1]
        lonsi = np.linspace(lons_m[0],lons[-1],nlons)
        latsi = np.linspace(lats_m[0],lats[-1],nlats)

        # also increase resolution of our edge lats/lons
        lats_e=regularbounds(latsi);
        lons_e=regularbounds(lonsi)
        lonsi, latsi = np.meshgrid(lonsi, latsi)
        # Smoothe data to increased resolution
        data = interp(data,lons,lats,lonsi,latsi)

    # Make edges into 2D meshed grid
#.........这里部分代码省略.........

def main():
    startTime = time.time()
    """Define the start and end date you want data extracted for:"""
    startYear=2009
    startMonth=10
    endYear=2012
    endMonth=12
    maxTries=3
    delay=10
    firstIteration=True
    lastIteration=False
    createFigure=False
    figureNumber=0
    USENETCDF4=True    # if false then use NETCDF3_CLASSIC


    """Name of output file to be created"""
    outputFile="NS8KM_obsSST_%s_to_%s.nc"%(startYear,endYear)
    if os.path.exists(outputFile): os.remove(outputFile)

    """Read the grid info from the grid file"""
    filename="/Users/trondkr/Projects/is4dvar/Grid/nordsjoen_8km_grid_hmax20m_v3.nc"
    mask_rho, lon_rho,lat_rho,grid_h = getGrid(filename)

    """Calculate the x,y grid coordinates"""
    (Mp,Lp)=lon_rho.shape
    X=np.arange(0,Mp,1)
    Y=np.arange(0,Lp,1)

    roms_Xgrid,roms_Ygrid=np.meshgrid(Y,X)

    """CoRTAD time is days since 1980/12/31 12:00:00"""
    mytime=getCortad.getCORTADtime()
    refDate=datetime.datetime(1981,12,31,12,0,0)

    """Have to convert the day of observation to the relative time used by ROMS
    which is 1948/1/1:00:00:00"""
    refDateROMS=datetime.datetime(1948,1,1,0,0,0)
    delta=refDate-refDateROMS
    daysSince1948to1980=delta.days

    """Find the start and end indexes to extract"""
    foundStart=False; foundEnd=False; startIndex=-9; endIndex=-9
    for index in xrange(len(mytime)):
        currentDate = refDateROMS + datetime.timedelta(days=float(mytime[index])+daysSince1948to1980)
        if foundStart is False:
            if currentDate.year==startYear:
                if currentDate.month==startMonth:
                    foundStart=True
                    startIndex=index
                    print "\n-----------------------------------------------"
                    print "Start date %s at index %s"%(currentDate,startIndex)

        if foundEnd is False:
            if currentDate.year==endYear:
                if currentDate.month==endMonth:
                    foundEnd=True
                    endIndex=index
                    print "FIXME : HARDCODING LAST INDEX !!!!!!!!!!!!!!!!!\n\n\n"
                    endIndex=1616
                    currentDate = refDateROMS + datetime.timedelta(days=float(mytime[endIndex])+daysSince1948to1980)
                    print "FIXME : HARDCODING LAST INDEX !!!!!!!!!!!!!!!!!\n\n\n"

                    print "End date %s at index %s"%(currentDate,endIndex)
    times=[i for i in range(startIndex,endIndex,1)]
    print "Created array of %s time-steps to iterate and extract data from"%(len(times))
    print "-----------------------------------------------\n"

    """Get the lomgitude-latitudes of the combination of tiles"""
    longitude, latitude, lonSST, latSST, indexes = getCortad.extractCoRTADLongLat(maxTries,
                                                                        delay,
                                                                        lon_rho.min(),
                                                                        lon_rho.max(),
                                                                        lat_rho.min(),
                                                                        lat_rho.max())
    indexes=np.asarray(indexes,dtype=np.int32)
    latitude  = np.flipud(latitude[indexes[3]:indexes[2]])
    longitude = longitude[indexes[0]:indexes[1]]

    """Loop over all times and store to file or make map"""
    polygon_data = getPolygon(lonSST[indexes[3]:indexes[2],indexes[0]:indexes[1]],
                              latSST[indexes[3]:indexes[2],indexes[0]:indexes[1]],
                              lon_rho,lat_rho)
    survey_time=[]

    for t in xrange(len(times)):
        
        """Open the files and check that NOAA is online"""
        cdf = getCortad.openCoRTAD(maxTries,delay)
        currentDate=refDateROMS + datetime.timedelta(days=int(mytime[times[t]])+daysSince1948to1980)

        
        """ Get the data for the current time"""
        filledSST = getCortad.extractCORTADSST("North Sea",times[t],cdf,indexes)

        """Interpolate the original values to the grid. This is the data that will be saved to file"""

        SSTi = mp.interp(np.flipud(filledSST),longitude,latitude,
                             lon_rho,lat_rho,checkbounds=False,masked=True,order=1)

#.........这里部分代码省略.........