def test_read_binary_subset_int(self):
        # INT16 = H
        fname = tempfile.mktemp()
        f = open(fname, 'w')
        ref = (self.x*10).astype('int16')
        f.write(ref)
        f.close()

        D = Data(None, None)
        f = open(fname, 'r')
        ny, nx = self.x.shape
        nt = 1

        # test 1: read entire file
        file_content = D._read_binary_subset2D(f, 2, ny=ny, nx=nx, xbeg=0, xend=nx, ybeg=0, yend=ny)
        d = np.reshape(np.asarray(struct.unpack('H'*ny*nx*nt, file_content)), (ny, nx))
        self.assertTrue(np.all(d-ref == 0.))

        # test 2: read subset with 1-values only
        ny1 = self.ymax - self.ymin
        nx1 = self.xmax - self.xmin
        nt1 = 1

        file_content = D._read_binary_subset2D(f, 2, ny=ny, nx=nx, xbeg=self.xmin, xend=self.xmax, ybeg=self.ymin, yend=self.ymax)
        d1 = np.reshape(np.asarray(struct.unpack('H'*ny1*nx1*nt1, file_content)), (ny1, nx1))
        self.assertTrue(np.all(d1 - ref[self.ymin:self.ymax, self.xmin:self.xmax] == 0.))

def test_mean_model():

    #The following code provides a routine that allows to validate the MeanModel() class
    print ('Jetzt gehts los')
    # generate some sample data ---
    x = Data(None, None)
    x.data = np.random.random((10,20,30))
    x.label='nothing'

    y = x.mulc(0.3)
    z = x.mulc(0.5)
    m = x.add(y).add(z).divc(3.)
    r = m.div(x)  # gives 0.6 as reference solution

    # generate Model instances and store Data objects as 'variables' ---
    dic_variables = ['var1', 'var2']
    X = Model(None, dic_variables, name='x', intervals='season')
    X.variables = {'var1': x, 'var2': x}
    Y = Model(None, dic_variables, name='y', intervals='season')
    Y.variables = {'var1': y, 'var2': y}
    Z = Model(None, dic_variables, name='z', intervals='season')
    Z.variables={'var1': z, 'var2': z}

    #... now try multimodel ensemble
    M=MeanModel(dic_variables,intervals='season')
    M.add_member(X)
    M.add_member(Y)
    M.add_member(Z)
    M.ensmean()  # calculate ensemble mean
    # print M.variables['var2'].div(x).data #should give 0.6
    npt.assert_equal(np.all(np.abs(1. - M.variables['var2'].div(x).data/0.6) < 0.00000001), True)

def get_T63_landseamask(shift_lon, mask_antarctica=True, area='land'):
    """
    get JSBACH T63 land sea mask
    the LS mask is read from the JSBACH init file

    area : str
        ['land','ocean']: When 'land', then the mask returned
        is True on land pixels, for ocean it is vice versa.
        In any other case, you get a valid field everywhere (globally)

    mask_antarctica : bool
        if True, then the mask is FALSE over Antarctica (<60S)
    """
    ls_file = get_data_pool_directory() \
        + 'variables/land/land_sea_mask/jsbach_T63_GR15_4tiles_1992.nc'
    ls_mask = Data(ls_file, 'slm', read=True, label='T63 land-sea mask',
                   lat_name='lat', lon_name='lon', shift_lon=shift_lon)
    if area == 'land':
        msk = ls_mask.data > 0.
    elif area == 'ocean':
        msk = ls_mask.data == 0.
    else:
        msk = np.ones(ls_mask.data.shape).astype('bool')

    ls_mask.data[~msk] = 0.
    ls_mask.data[msk] = 1.
    ls_mask.data = ls_mask.data.astype('bool')
    if mask_antarctica:
        ls_mask.data[ls_mask.lat < -60.] = False

    return ls_mask

 def setUp(self):
     D = Data(None, None)
     D.data = np.random.random((10, 20))
     lon = np.arange(-10.,10.)  # -10 ... 9
     lat = np.arange(-60., 50., 2.)  # -60 ... 48
     D.lon, D.lat = np.meshgrid(lon, lat)
     self.x = D

 def setUp(self):
     D = Data(None, None)
     tmp = np.random.random((55, 20))
     D.data = np.ma.array(tmp, mask=tmp!=tmp)
     lon = np.arange(-10.,10.)  # -10 ... 9
     lat = np.arange(-60., 50., 2.)  # -60 ... 48
     LON, LAT = np.meshgrid(lon, lat)
     D.lon = np.ma.array(LON, mask=LON!=LON)
     D.lat = np.ma.array(LAT, mask=LAT!=LAT)
     self.x = D

    def test_read_full_binary_file_double(self):
        # write binary test data
        fname = tempfile.mktemp()
        f = open(fname, 'w')
        f.write(self.x)
        f.close()

        D = Data(None, None)
        D.filename = fname
        ny, nx = self.x.shape
        D._read_binary_file(ny=ny, nx=nx, nt=1, dtype='double')
        self.assertTrue(np.all(D.data-self.x == 0.))

    def setUp(self):
        n=1000  # slows down significantly! constraint is percentile  test
        x = sc.randn(n)*100.  # generate dummy data
        self.D = Data(None, None)
        d=np.ones((n, 1, 1))
        self.D.data = d
        self.D.data[:,0,0]=x
        self.D.data = np.ma.array(self.D.data, mask=self.D.data != self.D.data)
        self.D.verbose = True
        self.D.unit = 'myunit'
        self.D.label = 'testlabel'
        self.D.filename = 'testinputfilename.nc'
        self.D.varname = 'testvarname'
        self.D.long_name = 'This is the longname'
        self.D.time = np.arange(n) + pl.datestr2num('2001-01-01')
        self.D.time_str = "days since 0001-01-01 00:00:00"
        self.D.calendar = 'gregorian'
        self.D.oldtime=False

        # generate dummy Model object
        data_dir = './test/'
        varmethods = {'albedo':'get_albedo()', 'sis': 'get_sis()'}
        self.model = models.Model(data_dir, varmethods, name='testmodel', intervals='monthly')

        sis = self.D.copy()
        sis.mulc(5., copy=False)
        sis.label='sisdummy'

        alb = self.D.copy()
        alb.label='albedodummy'

        # add some dummy data variable
        self.model.variables = {'albedo':alb, 'sis':sis}

    def get_temperature_2m(self, interval=None):
        """
        return data object of
        a) seasonal means for air temperature
        b) global mean timeseries for TAS at original temporal resolution
        """
        print 'Needs revision to support CMIP RAWDATA!!'
        assert False

        if interval != 'season':
            raise ValueError('Other data than seasonal not supported at the moment for CMIP5 data and temperature!')

        #original data
        filename1 = self.data_dir + 'tas/' + self.model + '/' + 'tas_Amon_' + self.model + '_' + self.experiment + '_ensmean.nc'

        force_calc = False

        if self.start_time is None:
            raise ValueError('Start time needs to be specified')
        if self.stop_time is None:
            raise ValueError('Stop time needs to be specified')

        s_start_time = str(self.start_time)[0:10]
        s_stop_time = str(self.stop_time)[0:10]

        tmp = pyCDO(filename1, s_start_time, s_stop_time, force=force_calc).seldate()
        tmp1 = pyCDO(tmp, s_start_time, s_stop_time).seasmean()
        filename = pyCDO(tmp1, s_start_time, s_stop_time).yseasmean()

        if not os.path.exists(filename):
            print 'WARNING: Temperature file not found: ', filename
            return None

        tas = Data(filename, 'tas', read=True, label=self._unique_name, unit='K', lat_name='lat', lon_name='lon', shift_lon=False)

        tasall = Data(filename1, 'tas', read=True, label=self._unique_name, unit='K', lat_name='lat', lon_name='lon', shift_lon=False)
        if tasall.time_cycle != 12:
            raise ValueError('Timecycle of 12 expected here!')

        tasmean = tasall.fldmean()
        retval = (tasall.time, tasmean, tasall)
        del tasall

        tas.data = np.ma.array(tas.data, mask=tas.data < 0.)

        return tas, retval

    def test_rasterize_init(self):
        x = Data(None, None)
        x._init_sample_object(ny=1, nx=272)
        x.lon = np.random.random(272)*10. + 5.  # 5 ... 15
        x.lat = np.random.random(272)*20. + 0.  # 0 ... 20

        lon = np.random.random((10,20))
        lat = np.random.random((30,20))

        with self.assertRaises(ValueError):
            x._rasterize(lon, lat, radius=0.1)
        lon = np.random.random((10,20))
        lat = np.random.random((10,20))

        with self.assertRaises(ValueError):
            x._rasterize(lon, lat, radius=None)

    def __init__(self, filename, gridfile, varname, read=False, **kwargs):
        """
        Parameters
        ----------

        filename : str
            filename of data file

        gridfile : str
            filename of grid definition file

        varname : str
            name of variable to handle

        read : bool
            specify if data should be read immediately
        """
        Data.__init__(self, filename, varname, **kwargs)
        self.gridfile = gridfile
        self.gridtype = 'unstructured'

    def setUp(self):
        #init Data object for testing
        n=4 #slows down significantly! constraint is percentile  test
        x = sc.randn(n)*100. #generate dummy data
        self.D = Data(None,None)
        d=np.ones((n,1,2))
        self.D.data = d
        self.D.data[:,0,0]=x
        self.D.data = np.ma.array(self.D.data,mask=self.D.data != self.D.data)
        self.D.verbose = True
        self.D.unit = 'myunit'

        self.D.time = np.arange(n) + pl.datestr2num('2001-01-01') - 1

    def test_read_binary_subset_Data_int(self):
        # binary data from subset in Data object

        # write binary test data
        fname = tempfile.mktemp()
        f = open(fname, 'w')

        tmp = (np.random.random(self.x.shape)*100.).astype('int16')
        f.write(tmp)
        f.close()

        D = Data(None, None)
        D.filename = fname
        ny, nx = self.x.shape

        latmin = self.lat[self.ymin]
        latmax = self.lat[self.ymax]
        lonmin = self.lon[self.xmin]
        lonmax = self.lon[self.xmax]

        D._read_binary_file(nt=1, dtype='int16', latmin=latmin, latmax=latmax, lonmin=lonmin, lonmax=lonmax, lat=self.lat, lon=self.lon)
        self.assertTrue(np.all(D.data-tmp[self.ymin:self.ymax+1,self.xmin:self.xmax+1] == 0.))

    def test_rasterize_data(self):
        """
        testdataset

        +---+---+---+
        |1.2|2.3|   |
        +---+---+---+
        |   |   |0.7|
        +---+---+---+
        |   |5.2|   |
        +---+---+---+
        """
        x = Data(None, None)
        x._init_sample_object(ny=1, nx=272)

        x.lon = np.asarray([2.25, 2.45, 1.8, 3.6])
        x.lat = np.asarray([11.9, 10.1, 10.2, 11.3])
        x.data = np.asarray([5.2, 2.3, 1.2, 0.7])

        # target grid
        lon = np.asarray([1.5, 2.5, 3.5])
        lat = np.asarray([10., 11., 12.])
        LON, LAT = np.meshgrid(lon, lat)

        # rasterize data

        # no valid data
        res = x._rasterize(LON, LAT, radius=0.000001, return_object=True)
        self.assertEqual(res.data.mask.sum(), np.prod(LON.shape))

        with self.assertRaises(ValueError):
            res = x._rasterize(LON, LAT, radius=0.000001, return_object=False)

        # check valid results
        res = x._rasterize(LON, LAT, radius=0.5, return_object=True)
        self.assertEqual(res.data[0,0], 1.2)
        self.assertEqual(res.data[0,1], 2.3)
        self.assertEqual(res.data[1,2], 0.7)
        self.assertEqual(res.ny*res.nx - res.data.mask.sum(), 4)

    def setUp(self):
        self.D = Data(None, None)
        self.D._init_sample_object(nt=1000, ny=1, nx=1)

        # generate dummy Model object
        data_dir = './test/'
        varmethods = {'albedo':'get_albedo()', 'sis': 'get_sis()'}
        self.model = models.Model(data_dir, varmethods, name='testmodel', intervals='monthly')

        sis = self.D.copy()
        sis.mulc(5., copy=False)
        sis.label='sisdummy'

        alb = self.D.copy()
        alb.label='albedodummy'

        # add some dummy data variable
        self.model.variables = {'albedo':alb, 'sis':sis}

 def setUp(self):
     # init Data object for testing
     n=100  # slows down significantly! constraint is percentile  test
     x = sc.randn(n)*100.  # generate dummy data
     self.D = Data(None, None)
     d=np.ones((n, 1, 1))
     self.D.data = d
     self.D.data[:,0,0]=x
     self.D.data = np.ma.array(self.D.data, mask=self.D.data != self.D.data)
     self.D.verbose = True
     self.D.unit = 'myunit'
     self.D.label = 'testlabel'
     self.D.filename = 'testinputfilename.nc'
     self.D.varname = 'testvarname'
     self.D.long_name = 'This is the longname'
     self.D.time = np.arange(n) + pl.datestr2num('2001-01-01') - 1
     self.D.time_str = "days since 0001-01-01 00:00:00"
     self.D.calendar = 'gregorian'
     self.D.cell_area = np.ones_like(self.D.data[0,:,:])

    def setUp(self):
        self.nx = 20
        self.ny = 10
        self.tempfile = tempfile.mktemp(suffix='.nc')
        self.gfile1 = tempfile.mktemp(suffix='.nc')
        self.gfile2 = tempfile.mktemp(suffix='.nc')
        self.gfile3 = tempfile.mktemp(suffix='.nc')
        self.x = Data(None, None)
        self.x._init_sample_object(nt=10, ny=self.ny, nx=self.nx)
        self.x.save(self.tempfile, varname='myvar')

        # generate some arbitrary geometry file
        F = NetCDFHandler()
        F.open_file(self.gfile1, 'w')
        F.create_dimension('ny', size=self.ny)
        F.create_dimension('nx', size=self.nx)
        F.create_variable('lat', 'd', ('ny', 'nx'))
        F.create_variable('lon', 'd', ('ny', 'nx'))
        F.assign_value('lat', np.ones((self.ny,self.nx)) * 5.)
        F.assign_value('lon', np.ones((self.ny,self.nx)) * 3.)
        F.close()

        F = NetCDFHandler()
        F.open_file(self.gfile2, 'w')
        F.create_dimension('ny', size=self.ny)
        F.create_dimension('nx', size=self.nx)
        F.create_variable('latitude', 'd', ('ny', 'nx'))
        F.create_variable('longitude', 'd', ('ny', 'nx'))
        F.assign_value('latitude', np.ones((self.ny,self.nx)) * 7.)
        F.assign_value('longitude', np.ones((self.ny,self.nx)) * 8.)
        F.close()

        F = NetCDFHandler()
        F.open_file(self.gfile3, 'w')
        F.create_dimension('ny', size=self.ny*2)
        F.create_dimension('nx', size=self.nx*3)
        F.create_variable('latitude', 'd', ('ny', 'nx'))
        F.create_variable('longitude', 'd', ('ny', 'nx'))
        F.assign_value('latitude', np.ones((self.ny*2,self.nx*3)) * 7.)
        F.assign_value('longitude', np.ones((self.ny*2,self.nx*3)) * 8.)
        F.close()

"""
This file is part of pyCMBS.
(c) 2012- Alexander Loew
For COPYING and LICENSE details, please refer to the LICENSE file
"""

"""
development script for pattern correlation analysis
"""

from pycmbs.diagnostic import PatternCorrelation
from pycmbs.data import Data
import numpy as np

import matplotlib.pyplot as plt

plt.close('all')
fname = '../pycmbs/examples/example_data/air.mon.mean.nc'

# generate two datasets
x = Data(fname, 'air', read=True)
xc = x.get_climatology(return_object=True)
yc = xc.copy()
yc.data = yc.data * np.random.random(yc.shape)*10.

PC = PatternCorrelation(xc, yc)
PC.plot()

plt.show()

class TestPycmbsBenchmarkingModels(unittest.TestCase):

    def setUp(self):
        self.D = Data(None, None)
        self.D._init_sample_object(nt=1000, ny=1, nx=1)

        # generate dummy Model object
        data_dir = './test/'
        varmethods = {'albedo':'get_albedo()', 'sis': 'get_sis()'}
        self.model = models.Model(data_dir, varmethods, name='testmodel', intervals='monthly')

        sis = self.D.copy()
        sis.mulc(5., copy=False)
        sis.label='sisdummy'

        alb = self.D.copy()
        alb.label='albedodummy'

        # add some dummy data variable
        self.model.variables = {'albedo':alb, 'sis':sis}

    def test_save_prefix_missing(self):
        m = self.model
        odir = tempfile.mkdtemp() + os.sep
        with self.assertRaises(ValueError):
            m.save(odir)

    def test_save_create_odir(self):
        m = self.model
        odir = tempfile.mkdtemp() + os.sep
        if os.path.exists(odir):
            os.system('rm -rf ' + odir)
        m.save(odir, prefix='test')
        self.assertTrue(os.path.exists(odir))
        os.system('rm -rf ' + odir)

    def test_save(self):
        m = self.model
        odir = tempfile.mkdtemp() + os.sep

        sisfile = odir + 'testoutput_SIS.nc'
        albfile = odir + 'testoutput_ALBEDO.nc'
        if os.path.exists(sisfile):
            os.remove(sisfile)
        if os.path.exists(albfile):
            os.remove(albfile)

        m.save(odir, prefix='testoutput')
        self.assertTrue(os.path.exists(sisfile))
        self.assertTrue(os.path.exists(albfile))

        if os.path.exists(sisfile):
            os.remove(sisfile)
        if os.path.exists(albfile):
            os.remove(albfile)
        os.system('rm -rf ' + odir)

    def test_cmip5_init_singlemember(self):
        data_dir = tempfile.mkdtemp()

        # invalid model identifier
        with self.assertRaises(ValueError):
            M = models.CMIP5RAW_SINGLE(data_dir, 'MPI-M:MPI-ESM-LR1', 'amip', {}, intervals='monthly')
        with self.assertRaises(ValueError):
            M = models.CMIP5RAW_SINGLE(data_dir, 'MPI-M:MPI-ESM-LR#1#2', 'amip', {}, intervals='monthly')
        M1 = models.CMIP5RAW_SINGLE(data_dir, 'MPI-M:MPI-ESM-LR#1', 'amip', {}, intervals='monthly')
        M2 = models.CMIP5RAW_SINGLE(data_dir, 'MPI-M:MPI-ESM-LR#728', 'amip', {}, intervals='monthly')
        self.assertEqual(M1.ens_member, 1)
        self.assertEqual(M2.ens_member, 728)

    def test_cmip5_singlemember_filename(self):
        data_dir = tempfile.mkdtemp()

        # generate testfile
        testfile = data_dir + os.sep + 'MPI-M' + os.sep + 'MPI-ESM-LR' + os.sep + 'amip' + os.sep + 'mon' + os.sep + 'atmos' + os.sep + 'Amon' + os.sep + 'r1i1p1' + os.sep + 'ta' + os.sep + 'ta_Amon_MPI-ESM-LR_amip_r1i1p1_197901-200812.nc'
        os.makedirs(os.path.dirname(testfile))
        os.system('touch ' + testfile)
        self.assertTrue(os.path.exists(testfile))

        M = models.CMIP5RAW_SINGLE(data_dir, 'MPI-M:MPI-ESM-LR#1', 'amip', {}, intervals='monthly')
        f = M.get_single_ensemble_file('ta', mip='Amon', realm='atmos')
        self.assertTrue(os.path.exists(f))
        self.assertEqual(f, testfile)

def xxxxtest_median_model():
    x = Data(None, None)
    x.label = 'nothing'
    d = np.random.random((100, 1, 1))
    x.data = np.ma.array(d, mask= d!=d)

    # odd number and no masked values
    a = x.copy()
    a.data[:, 0, 0] = 1.
    b = x.copy()
    b.data[:, 0, 0] = 3.
    c = x.copy()
    c.data[:, 0, 0] = 2.
    d = x.copy()
    d.data[:, 0, 0] = 5.
    e = x.copy()
    e.data[:, 0, 0] = 4.

    m = MedianModel()
    m.add_member(a)
    m.add_member(b)
    m.add_member(c)
    m.add_member(d)
    m.add_member(e)
    m.ensmedian()

    # should give the value of 3. for all timesteps

    del m

    # even number and no masked values
    a = x.copy()
    a.data[:, 0, 0] = 1.
    b = x.copy()
    b.data[:, 0, 0] = 3.
    c = x.copy()
    c.data[:, 0, 0] = 2.
    d = x.copy()
    c.data[:, 0, 0] = 4.

    m = MedianModel()
    m.add_member(a)
    m.add_member(b)
    m.add_member(c)
    m.add_member(d)

    m.ensmedian()

    # should give the value of 2.5 for all timesteps
    del m

class TestEOF(TestCase):

    def setUp(self):
        #init Data object for testing
        n=4 #slows down significantly! constraint is percentile  test
        x = sc.randn(n)*100. #generate dummy data
        self.D = Data(None,None)
        d=np.ones((n,1,2))
        self.D.data = d
        self.D.data[:,0,0]=x
        self.D.data = np.ma.array(self.D.data,mask=self.D.data != self.D.data)
        self.D.verbose = True
        self.D.unit = 'myunit'

        self.D.time = np.arange(n) + pl.datestr2num('2001-01-01') - 1


    def test_eof_analysis(self):
        #test of EOF
        #example taken from:
        #   http://www.atmos.washington.edu/~dennis/552_Notes_4.pdf , page 80

        #assign sampled data
        D = self.D.copy()
#        d1 = np.array([2,4,-6,8])
#        d2 = np.array([1,2,-3,4])
        #d1 = np.array([2,1])
        #d2 = np.array([4,2])
        #d3 = np.array([-6,-3])
        #d4 = np.array([8,4])

        #D.data[:,0,0] = d1; D.data[:,0,1] = d2
        #D.data[:,0,2] = d3; D.data[:,0,3] = d4

#        D.data[:,0,0] = d1
#        D.data[:,0,1] = d2
#
#        E = EOF(D,cov_norm=False) #do not normalize covariance matrix, as this is also not done in example
#        print E.C #covariance matrix


#        shape of EOF is wrong !!!!!!!
#
#        something is not really working here !!!!
#
#
#        irgendwie ist hier ein problem
#        warum einmal 4x4 und einmal 2x2 ???

#        print 'eigval'
#        print E.eigval
#
#        print 'eigvec'
#        print E.eigvec






    pass