.. _tutorial_vima_scidata2: ============================== How to handle IMA science data ============================== Due to highly complex data structure of the IMA data (too many modes, time series of 4-dimensional, many processing levels, different operations, similar sensors on MEX & VEX but not fully compatible operations and processing, ...) the access to the IMA science data is highly complicated. The ``irfpy.aspera`` library give a quick access to such dataset. It is noted that the "final form of the matrix" depends highly on the users. Thus, the interface of the data access should be as simple and flexible as possible. Getting the matrix at the specific time ======================================= To get the full matrix (192s data), you may just use the following syntax. >>> import datetime >>> from irfpy.vima import scidata_util # Use from irfpy.mima import scidata_util for MEX. >>> tstart, data = scidata_util.get3d(datetime.datetime(2008, 1, 9, 14, 30)) The returned ``tstart`` is the time of the start of the data. The ``data`` includes the VEX/IMA data. The ``data.matrix`` is the matrix that you want. >>> print(data.matrix) This is (32, 16, 96, 8) shaped array. Getting the time series matrix ============================== To get the time series of the matrix, you may use the following syntax. >>> import datetime >>> from irfpy.vima import scidata_util >>> timeser = scidata_util.getarray2d(datetime.datetime(2008, 1, 9, 14, 30), (datetime.datetime(2008, 1, 9, 15, 30))) >>> print(timeser) A bit more in details ===================== ASPERA/IMA data structure ------------------------- The ASPERA/IMA sensor can measure ions with mass, energy and direction discrimilated. At a certain time, IMA measures ions at a specific energy and polar angle. Simultaneously, the mass and azimuth angle can be discrimimated. It means that the data in IMA at a certain time has (M=32, A=16) shaped array. The accumulation time is 125 ms. Then, the energy is swept. Usually, energy has 96 steps with 12 sec. Considering this, the 12-sec resolution data with (M=32, A=16, E=96) can be obtained. The energy sweep is repeated with different polar angle scan. 16 steps (=192 s) are sweeped. The resulting dataset is the 192-sec resolution data with (M=32, A=16, E=96, P=16). This 192-sec resolution data is transferred as is if the mode is 24. In this mode, the data is 192-sec resolution data with (M=32, A=16, E=96, P=16). For MEX, mode=24 is used near the pericenter. If the mode is 25, the neighboring polar angle is added up on board, and sent to the ground. In this mode, the data is 192-sec resolution data with (M=32, A=16, E=96, P=8). For VEX, mode=25 has been used for the beginning of the mission. Later the mode=24 is used, but with 1 count reduction, for VEX. In summary, the data at a certain time can be considered as (depending on the way of data analysis) 1. 3D (azimuth, energy, polar + mass) data with time resolution of 192s 2. 2D (azimuth, energy + mass) data with time resolution of 12s (mode 24) 3. 2D (azimuth, energy + mass) data with time resolution of 24s (mode 25) 4. 1D (azimuth) data with resolution of 125 ms at specific energy and (mode 24) 5. 1D (azimuth) data with resolution of 250 ms at specific energy and (mode 25) Each of them is indeed a time series. More advanced dataset --------------------- For IMA, there are other "operation", such as 32 second mode, and fast-sweep mode (24-s resolution). 32 second mode has not used so that IRFPY will not support this. Fast sweep mode was prepared for Phobos flyby for MEX. This is intended to be implemented into IRPFY. Object oriented dataset ----------------------- For ``irfpy.aspera``, each of these data is represented by a specific class in an object oriented way. Below we review the class used for each dataset. =================== ============================================== ======================================================== Data At specific time Time series =================== ============================================== ======================================================== 3D (Az, En, Po; Ma) :class:`irfpy.imacommon.imascipac.CntMatrix` :class:`irfpy.imacommon.imascipac.TimeSeriesCntMatrix` 2D (Az, En; Ma) :class:`irfpy.imacommon.imascipac.CntMatrix2D` :class:`irfpy.imacommon.imascipac.TimeSeriesCntMatrix2D` 1D (Az; Ma) N/A (To be developed) N/A (To be developed) =================== ============================================== ======================================================== 3D (azimuth, energy, polar + mass) .................................. To represent the 3D distribution data at a certain time, :class:`irfpy.imacommon.imascipac.CntMatrix` class is used. The object contains the count rate matrix, mode, and observation time. Count rate matrix is a masked array with the shape corresponding to the mode (for Mode 24, the shape is (32, 16, 96, 16)). 2D (azimuth, energy + mass) data ................................ To represent the 2D distribution at a ceratin time, :class:`irfpy.imacommon.imascipac.CntMatrix2D` class is used. The object contains the count rate matrix, mode, polar index, and observation time. Count rate matrix is a masked array with the shape corresponding to the mode (but both mode 24 and 25, the shape is (32, 16, 96)). The polar index is a pair of (start, stop), both inclusive. Time series of the data ....................... Time series of 3D data is :class:`irfpy.imacommon.imascipac.TimeSeriesCntMatrix`. Use :class:`irfpy.vima.scidata_util.getarray3d` to get the time series data. Time series of 2D data is :class:`irfpy.imacommon.imascipac.TimeSeriesCntMatrix2D`. Use :class:`irfpy.vima.scidata_util.getarray2d` to get the time series data. To read (VEX/IMA) ----------------- The first step to read the VEX/IMA is getting time series of 2D data. Reading time series of 2D data .............................. Use :func:`irfpy.vima.scidata_util.getarray2d` function. With time intervals as argument, the function will return the time series of 2D data (:class:`irfpy.imacommon.imascipac.TimeSeriesCntMatrix2D`. >>> import datetime >>> from irfpy.vima import scidata_util >>> t0 = datetime.datetime(2012, 10, 5, 0, 1, 23) >>> t1 = datetime.datetime(2012, 10, 6, 4, 59, 43) >>> dataset = scidata_util.getarray2d(t0, t1) The dataset represent the time series of VEX/IMA 2D data. Let's look at the `dataset`. >>> print(dataset) You may find that there are 3362 2-D data. You can use the :meth:`getobstime() ` method to get the list of observation time. >>> obstime = dataset.getobstime() >>> print(obstime) [datetime.datetime(2012, 10, 5, 0, 1, 11, 198276), datetime.datetime(2012, 10, 5, 0, 1, 23, 198281), datetime.datetime(2012, 10, 5, 0, 1, 36, 167030), datetime.datetime(2012, 10, 5, 0, 1, 48, 167035), datetime.datetime(2012, 10, 5, 0, 2, 0, 167040), ...] You can get the array expression of the dataset by :meth:`getarray() ` method. >>> dataarray = dataset.getarray() >>> print(dataarray.shape) (32, 16, 96, 3362) Note on the start time '''''''''''''''''''''' It could be slightly minor issue, but the first data in the example above is slightly earlier than the given time. Given time is ``t0=2012-10-05T00:01:23``, but the first data is observed by ``obstime[0]=2012-10-05T00:01:11``. This is because the data at ``t0`` might be included in the data at index 0. This could produce a problem if ``t0`` falls at a large data gap. For example, >>> tg0 = datetime.datetime(2006, 12, 5, 0, 0, 0) >>> tg1 = datetime.datetime(2006, 12, 6, 0, 0, 0) >>> dataset2 = scidata_util.getarray2d(tg0, tg1) >>> obstime2 = dataset2.getobstime() >>> print(obstime2) [datetime.datetime(2006, 12, 4, 9, 41, 0, 48676), datetime.datetime(2006, 12, 5, 5, 21, 31, 41914), datetime.datetime(2006, 12, 5, 5, 21, 55, 41914), datetime.datetime(2006, 12, 5, 5, 22, 19, 41925), datetime.datetime(2006, 12, 5, 5, 22, 43, 41926), ...] As you see, the first data is observed 14 hour before the given time. Thus, this is users' responsibility to check if the first data (index=0) is valid for the data analysis or not. One idea is to use the following syntax. >>> if t0 - obstime2[0] > datetime.timedelta(seconds=12): ... dataset2 = dataset2.clipped(obstime2[1], obstime2[-1]) >>> obstime2 = dataset2.getobstime() >>> print(obstime2) [datetime.datetime(2006, 12, 5, 5, 21, 31, 41914), datetime.datetime(2006, 12, 5, 5, 21, 55, 41914), datetime.datetime(2006, 12, 5, 5, 22, 19, 41925), datetime.datetime(2006, 12, 5, 5, 22, 43, 41926), ...] Reading 2D data at a certain time ................................. You can get the single 2-D data using the index of the time series data, iteration, and :meth:`get() ` method. In any cases, the obtained value is a tuple of (`datetime`, :class:`irfpy.imacommon.imascipac.CntMatrix2D`). * Using index >>> data0 = dataset[0] # Return the first data. >>> print(data0) (datetime.datetime(2012, 10, 5, 0, 1, 11, 198276), ) * Using iteration >>> for obstime, data in dataset: ... print(obstime, data) 2012-10-05 00:01:11.198276 <@2012-10-05T00:01:11.198276:MOD=24 >>24<<:POL=14/14:CNTmax=13> 2012-10-05 00:01:23.198281 <@2012-10-05T00:01:23.198281:MOD=24 >>24<<:POL=15/15:CNTmax=5> 2012-10-05 00:01:36.167030 <@2012-10-05T00:01:36.167030:MOD=24 >>24<<:POL=00/00:CNTmax=2> 2012-10-05 00:01:48.167035 <@2012-10-05T00:01:48.167035:MOD=24 >>24<<:POL=01/01:CNTmax=4> ... * Using ``get`` method. >>> data_at_0100 = dataset.get(datetime.datetime(2012, 10, 5, 1)) # Data at 1 oclock. >>> print(data_at_0100) (datetime.datetime(2012, 10, 5, 0, 59, 48, 230124), ) Now let's look at the :class:`irfpy.imacommon.imascipac.CntMatrix2D` object more carefully. >>> dat2d = data0[1] >>> print(dat2d) <@2012-10-05T00:01:11.198276:MOD=24:POL=14/14:CNTmax=13> You may see that the ``CntMatrix2D`` contains the information on the counts and some other misc data. The attributes below includes them. * t Time in datetime object * matrix Matrix of count rate. Shape depends on the mode, but for nominal modes it is (M32,A16,E96). * polar Tuple providing the start polar index and the stop polar index (inclusive) * mode Mode object (:class:`irfpy.imacommon.imamode.Mode`) See example below. >>> print(dat.t) 2012-10-05 00:01:11.198276 >>> print(dat.polar) [14, 14] >>> print(dat.mode) >>> print(dat.matrix) [[[0 0 0 ..., 0 0 0] [0 0 0 ..., 0 0 0] [0 0 0 ..., 0 0 0] ..., [0 0 0 ..., 0 0 0] [0 0 0 ..., 0 0 0] [0 0 0 ..., 0 0 0]] ... >>> print(dat.matrix.shape) (32, 16, 96) Time series of 3D data ...................... .. Time series of 2-D data can be converted to the time series of 3-D data. .. .. Just to repeat, we first get the time series of 2-D data as above: .. .. >>> import datetime .. >>> from irfpy.vima import scidata_util .. >>> t0 = datetime.datetime(2012, 10, 5, 0, 1, 23) .. >>> t1 = datetime.datetime(2012, 10, 6, 4, 59, 43) .. >>> dataset = scidata_util.getarray2d(t0, t1) .. >>> print(dataset) .. .. .. The ``dataset`` is the time series of 2-D data. .. Now the ``dataset`` is converted to the time series of 3-D data. .. .. >>> from irfpy.imacommon import imascipac .. >>> dataset3d = imascipac.convert2to3(dataset) Time series of 3-D data can be obtained by using :func:`irfpy.vima.scidata_util.getarray3d` function. >>> import datetime >>> from irfpy.vima import scidata_util >>> t0 = datetime.datetime(2012, 10, 5, 0, 1, 23) >>> t1 = datetime.datetime(2012, 10, 6, 4, 59, 43) >>> dataset3d = scidata_util.getarray3d(t0, t1) The ``dataset3d`` is an object of :class:`irfpy.imacommon.imascipac.TimeSeriesCntMatrix`, as seen below. >>> print(dataset3d) For the given time range, you can realize 211 3-D data is available. Of course the functionality works for the mode 25 data. >>> t0 = datetime.datetime(2008, 10, 5, 2, 1, 23) >>> t1 = datetime.datetime(2008, 10, 5, 8, 59, 43) .. >>> dataset = scidata_util.getarray2d(t0, t1) .. >>> dataset3d = imascipac.convert2to3(dataset) >>> dataset3d = scidata_util.getarray3d(t0, t1) >>> print(dataset3d) You can look at the observation time as below; the time resolution is 24 s for this case. >>> print(dataset3d.getobstime()) [datetime.datetime(2008, 10, 5, 0, 38, 38, 833861), datetime.datetime(2008, 10, 5, 6, 7, 3, 868441), datetime.datetime(2008, 10, 5, 6, 10, 15, 868476), datetime.datetime(2008, 10, 5, 6, 13, 27, 837259), datetime.datetime(2008, 10, 5, 6, 16, 39, 831038), datetime.datetime(2008, 10, 5, 6, 19, 51, 799822), datetime.datetime(2008, 10, 5, 6, 23, 3, 793610), ... You can get the array expression of the dataset by :meth:`getarray() ` method. >>> dataarray = dataset3d.getarray() >>> print(dataarray.shape) (32, 16, 96, 8, 54) Specific 3D data ................ To get the 3D data of the specific time, you can use the index, or :meth:`get` method. >>> print(dataset3d[15]) (datetime.datetime(2008, 10, 5, 6, 51, 51, 881401), ) >>> dat3d15 = dataset3d[15][1] # Index 0 is time, and index 1 is the data (CntMatrix object) >>> print(dat3d15.mode) The :attr:`matrix` contains the count rate data. .. code-block:: py >>> print(dat3d15.matrix) [[[[ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] ..., [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 1. 0. ..., 0. 0. 0.]] [[ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] ..., [ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] ... >>> print(dat3d15.matrix.shape) (32, 16, 96, 8) Processing data --------------- There are some functions that can convert the object-oriented dataset to a simple numpy array. :func:`irfpy.imacommon.imascipac.timeseries2d_to_energytime_simple` will convert the :class:`irfpy.imacommon.imascipac.TimeSeriesCntMatrix2D` to (96, N) numpy array. This is useful to make the E-t diagram. >>> import datetime >>> from irfpy.vima import scidata_util >>> t0 = datetime.datetime(2012, 10, 5, 0, 1, 23) >>> t1 = datetime.datetime(2012, 10, 6, 4, 59, 43) >>> dataset = scidata_util.getarray2d(t0, t1) >>> print(dataset) >>> from irfpy.imacommon import imascipac >>> etmatrix = imascipac.timeseries2d_to_energytime_simple(dataset) # etmatrix is numpy array. >>> print(etmatrix) [[ 0. 0. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.] ..., [ 0. 1. 0. ..., 0. 0. 0.] [ 0. 0. 0. ..., 0. 0. 0.]] >>> print(etmatrix.shape) (96, 3362)