mraspaud/fy3_mersi_reader.py

## fy3_mersi_reader.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-

# Copyright (c) 2015 Martin Raspaud

# Author(s):

#   Martin Raspaud <martin.raspaud@smhi.se>

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""Read MERSI data from framed raw data.
"""


from datetime import datetime, timedelta
import numpy as np
import sys
from reed import rs
from trollimage.image import Image
from bitstring import BitArray, BitStream
import logging

logger = None

#from zfec import Decoder
#from trollcast.formats.reedsolomon import *
#from reedsolo import RSCodec

table = np.array([0xff, 0x48, 0x0e, 0xc0, 0x9a, 0x0d, 0x70, 0xbc, 0x8e, 0x2c,
                  0x93, 0xad, 0xa7, 0xb7, 0x46, 0xce, 0x5a, 0x97, 0x7d, 0xcc,
                  0x32, 0xa2, 0xbf, 0x3e, 0x0a, 0x10, 0xf1, 0x88, 0x94, 0xcd,
                  0xea, 0xb1, 0xfe, 0x90, 0x1d, 0x81, 0x34, 0x1a, 0xe1, 0x79,
                  0x1c, 0x59, 0x27, 0x5b, 0x4f, 0x6e, 0x8d, 0x9c, 0xb5, 0x2e,
                  0xfb, 0x98, 0x65, 0x45, 0x7e, 0x7c, 0x14, 0x21, 0xe3, 0x11,
                  0x29, 0x9b, 0xd5, 0x63, 0xfd, 0x20, 0x3b, 0x02, 0x68, 0x35,
                  0xc2, 0xf2, 0x38, 0xb2, 0x4e, 0xb6, 0x9e, 0xdd, 0x1b, 0x39,
                  0x6a, 0x5d, 0xf7, 0x30, 0xca, 0x8a, 0xfc, 0xf8, 0x28, 0x43,
                  0xc6, 0x22, 0x53, 0x37, 0xaa, 0xc7, 0xfa, 0x40, 0x76, 0x04,
                  0xd0, 0x6b, 0x85, 0xe4, 0x71, 0x64, 0x9d, 0x6d, 0x3d, 0xba,
                  0x36, 0x72, 0xd4, 0xbb, 0xee, 0x61, 0x95, 0x15, 0xf9, 0xf0,
                  0x50, 0x87, 0x8c, 0x44, 0xa6, 0x6f, 0x55, 0x8f, 0xf4, 0x80,
                  0xec, 0x09, 0xa0, 0xd7, 0x0b, 0xc8, 0xe2, 0xc9, 0x3a, 0xda,
                  0x7b, 0x74, 0x6c, 0xe5, 0xa9, 0x77, 0xdc, 0xc3, 0x2a, 0x2b,
                  0xf3, 0xe0, 0xa1, 0x0f, 0x18, 0x89, 0x4c, 0xde, 0xab, 0x1f,
                  0xe9, 0x01, 0xd8, 0x13, 0x41, 0xae, 0x17, 0x91, 0xc5, 0x92,
                  0x75, 0xb4, 0xf6, 0xe8, 0xd9, 0xcb, 0x52, 0xef, 0xb9, 0x86,
                  0x54, 0x57, 0xe7, 0xc1, 0x42, 0x1e, 0x31, 0x12, 0x99, 0xbd,
                  0x56, 0x3f, 0xd2, 0x03, 0xb0, 0x26, 0x83, 0x5c, 0x2f, 0x23,
                  0x8b, 0x24, 0xeb, 0x69, 0xed, 0xd1, 0xb3, 0x96, 0xa5, 0xdf,
                  0x73, 0x0c, 0xa8, 0xaf, 0xcf, 0x82, 0x84, 0x3c, 0x62, 0x25,
                  0x33, 0x7a, 0xac, 0x7f, 0xa4, 0x07, 0x60, 0x4d, 0x06, 0xb8,
                  0x5e, 0x47, 0x16, 0x49, 0xd6, 0xd3, 0xdb, 0xa3, 0x67, 0x2d,
                  0x4b, 0xbe, 0xe6, 0x19, 0x51, 0x5f, 0x9f, 0x05, 0x08, 0x78,
                  0xc4, 0x4a, 0x66, 0xf5, 0x58], dtype=np.uint8)


def decode_str(arr):
    decoder = np.tile(table, np.ceil(len(arr) / 255.0))
    return (np.fromstring(arr, np.uint8) ^ decoder[:len(arr)]).tostring()


def show(data, filename=None, stripes=0):
    """Show the stetched data.
    """

    norm_arr = np.array((data - data.min()) * 255.0 /
                        (data.max() - data.min())).astype(np.uint8)
    if stripes:
        r_norm_arr = np.array(norm_arr)
        r_norm_arr[:, ::stripes] = np.uint8(255 - r_norm_arr[:, ::stripes])
        img = Image((r_norm_arr, norm_arr, norm_arr), mode="RGB")
        img.stretch("crude")
    else:
        from PIL import Image as pil
        img = pil.fromarray(norm_arr)

    #img = Image((ch14, ch12, ch09), mode="RGB")
    #img = Image((ch14[:, 38:], ch12[:, 38:], ch09[:, 38:]), mode="RGB")

    if filename:
        img.save(filename)
    else:
        img.show()


def dec12to16(data):
    arr12 = data.astype(np.uint16).flat
    new_shape = list(data.shape[:-1]) + [(data.shape[-1] * 8) / 12]
    arr16 = np.zeros(new_shape, dtype=np.uint16)
    arr16.flat[::2] = np.left_shift(arr12[::3], 4) + \
        np.right_shift((arr12[1::3]), 4)
    arr16.flat[1::2] = np.left_shift((arr12[1::3] & 15), 8) + arr12[2::3]
    return arr16


#############################################################

from datetime import datetime, timedelta
import numpy as np
import sys
#from reed import rs
#from zfec import Decoder
#from trollcast.formats.reedsolomon import *
#from reedsolo import RSCodec, ReedSolomonError


vcdu_type = np.dtype([('version', '>u2'),
                      ('vcdu count', '>u1', (3, )),
                      ('replay', '>u1'),
                      ('insert', '>u2'),
                      ('mpdu_hdr', '>u2')])

ccsds_fy_hdr = np.dtype([('ccsds_version', '>u2'),
                         ('sequence', '>u2'),
                         ('packet_length', '>u2'),
                         ('days', '>u2'),
                         ('milliseconds', '>u4')])


def cadu_info(line):
    print "=" * 3, "Cadu Info", "=" * 3
    print line["vcdu count"][0] * 256 * 256 + line["vcdu count"][1] * 256 + line["vcdu count"][2]
    print "vcid", line["version"] & (2 ** 6 - 1)
    print "spacecraft_id", (line["version"] >> 6) & (2 ** 8 - 1)
    print "version", (line["version"] >> 14) & 3
    if line["insert"]:
        print "Data crypted"
    else:
        print "Data open"
    print "offset", line["mpdu_hdr"] & (2 ** 11 - 1)


def read_vcdu_hdr(packets):
    for packet in packets:
        line = np.fromstring(packet, vcdu_type, count=1)[0], packet[-882:]
        # cadu_info(line[0])
        # raw_input()
        yield line


def reed_decode(packets):
    codec = rs.RSCoder(255, 223)
    res = bytearray("\0" * 892)
    #codec = RSCodec(32)

    for packet in packets:
        # do the decoding!
        yield packet[:892]
        # for i in range(4):
        #     #print i
        #     try:
        #         print len(packet[i::4])
        #         print codec.verify(packet[i::4])
        #         decoded = codec.decode(packet[i::4])
        #         if not decoded:
        #             raise ValueError
        #         #print decoded
        #         res[i::4] = decoded
        #         #print "decoded"
        #         #raw_input()
        #     except IOError:
        #         print "uncorrectable"
        #         res[i::4] = packet[i:892:4]
        #     except ReedSolomonError as err:
        #         #print "uncorrectable:", err
        #         res[i::4] = packet[i:892:4]
        # print str(res[:892]) == packet[:892]
        # yield str(res[:892])


def pn_decode(packets, dlen=0):
    for packet in packets:
        if dlen < 0:
            yield packet
        else:
            yield packet[:dlen] + decode_str(packet[dlen:])


def frame_sync(packets):
    """Phony frame_sync.
    """
    for packet in packets:
        yield packet[4:]


def file_reader(filename, plen=1024):
    with open(filename) as fd:
        while True:
            packet = fd.read(plen)
            if packet:
                yield packet
            else:
                return


def sort_vc(packets):
    channels = {}
    cnt = 0
    for hdr, packet in packets:
        vcid = hdr["version"] & (2 ** 6 - 1)
        channels.setdefault(vcid, []).append((hdr, packet))
        if cnt % 1000 == 0:
            logger.debug(str(cnt))
        cnt += 1
    return channels


def gen_strips(data):
    frame_sync = "\xff\xe0\x01\xff\xe0\x01\xff\xe0\x01"

    bla = BitStream(bytes=data[:55000])
    offset = bla.find(BitArray(bytes=frame_sync))[0] % 8
    logger.debug("Bit offset %d", offset)
    if offset != 0:
        array = np.fromstring(data, dtype=np.uint8)
        data = (((array[:-1] & (2**(8 - offset) - 1)) << offset) + (array[1:] >> (8 - offset))).astype(np.uint8).tostring()
    prev = 0
    while True:
        try:
            idx = data.index(frame_sync, prev + len(frame_sync))
        except ValueError:
            break
        yield data[prev:idx]
        prev = idx


def read5(filename, plen, pn_start):
    channels = sort_vc(
        read_vcdu_hdr(reed_decode(pn_decode(frame_sync(file_reader(filename, plen=plen)), pn_start))))
    logger.debug('Channel keys: %s', str(channels.keys()))

    return channels

if __name__ == '__main__':

    import argparse

    parser = argparse.ArgumentParser()
    parser.add_argument("-c", "--cadu-len", help="Length of the cadu package in bytes (usually 1024 bytes)",
                        type=int, default=1024)
    parser.add_argument("-p", "--pseudo-noise-start", help="Start of the pseudo-noise chunk, usually 0. Default is -1 for no decoding.",
                        type=int, default=-1)
    parser.add_argument("-v", "--verbose", help="print debug messages too",
                        action="store_true")
    parser.add_argument("filename", help="name of the file to read")
    opts = parser.parse_args()

    handler = logging.StreamHandler()
    handler.setFormatter(logging.Formatter("[%(levelname)s: %(asctime)s :"
                                           " %(name)s] %(message)s",
                                           '%Y-%m-%d %H:%M:%S'))
    if opts.verbose:
        loglevel = logging.DEBUG
    else:
        loglevel = logging.INFO
    handler.setLevel(loglevel)
    logging.getLogger('').setLevel(loglevel)
    logging.getLogger('').addHandler(handler)
    logger = logging.getLogger("fy3_mersi_reader")


    chs = read5(opts.filename, opts.cadu_len, opts.pseudo_noise_start)

    logger.debug(str(chs.keys()))

    data = "".join(dat[1] for dat in chs[3])

    res = []
    low_res = []
    truc = []
    rest = []
    truc_str = []
    logger.info('Stripping data')
    for strip in gen_strips(data):
        if len(strip) == 12345: #12345:
            a = np.fromstring(strip, dtype=np.uint8)
            res.append(dec12to16(a))

        elif len(strip) == 3129:
            a = np.fromstring(strip, dtype=np.uint8)
            low_res.append(dec12to16(a))
            truc.append(a)
        elif len(strip) == 54531:
            a = np.fromstring(strip[:3129], dtype=np.uint8)
            low_res.append(dec12to16(a))
            truc_str.append(strip[3129:])
            rest.append(dec12to16(np.fromstring(strip[3129:], dtype=np.uint8)))
    truc = np.vstack(truc)
    low_res = np.vstack(low_res)
    res = np.vstack(res)
    restarr = np.vstack(rest)

    #for i in range(100):
        #a = truc[150 * i, 12:30].astype(np.uint32)
        #a = np.take(a, np.arange(6) + 3, mode="wrap")
        # milli = ((a[3] & 15) << 28) + ((a[4] >> 4) << 24) + \
        #    ((a[1] & 15) << 20) + (a[2] << 12) + (a[0] << 4) + a[1] >> 4
        #milli = (((((a[3] >> 4) << 12) + a[3]) << 12) + a[2])

        #milli = 0
        #day = ((a[4] & 15) << 12) + (a[5] << 4) + (a[3] >> 4)

        # print a
        # half = (((a[:-1] & 15) << 4) + (a[1:] >> 4)).astype(np.uint8)
        # print half
        # print dec12to16(a)
        # day = (half[4] << 8) + half[5]
        # milli = 0

        #day = (a[0] << 8) + a[1]
        #milli = (a[2] << 24) + (a[3] << 16) + (a[4] << 8) + a[5]

        #thedate = datetime(
        #    2000, 1, 1) + timedelta(days=int(day) - 1, milliseconds=int(milli))
        #print thedate

    #print "waiting..."
    #raw_input()

    res_chn = ["01", "02", "03", "04", "05"]
    low_res_chn = ["01", "02", "03", "04", "05"]

    ch01 = res[res[:, 6] < 40, :]
    ch02 = res[(res[:, 6] >= 40) & (res[:, 6] < 80), :]
    ch03 = res[(res[:, 6] >= 80) & (res[:, 6] < 120), :]
    ch04 = res[(res[:, 6] >= 120) & (res[:, 6] < 160), :]
    ch05 = res[(res[:, 6] >= 160) & (res[:, 6] < 200), :]

    ch06 = low_res[low_res[:, 6] < 210, :]
    ch07 = low_res[(low_res[:, 6] >= 210) & (low_res[:, 6] < 220), :]
    ch08 = low_res[(low_res[:, 6] >= 220) & (low_res[:, 6] < 230), :]
    ch09 = low_res[(low_res[:, 6] >= 230) & (low_res[:, 6] < 240), :]
    ch10 = low_res[(low_res[:, 6] >= 240) & (low_res[:, 6] < 250), :]
    ch11 = low_res[(low_res[:, 6] >= 250) & (low_res[:, 6] < 260), :]
    ch12 = low_res[(low_res[:, 6] >= 260) & (low_res[:, 6] < 270), :]
    ch13 = low_res[(low_res[:, 6] >= 270) & (low_res[:, 6] < 280), :]
    ch14 = low_res[(low_res[:, 6] >= 280) & (low_res[:, 6] < 290), :]
    ch15 = low_res[(low_res[:, 6] >= 290) & (low_res[:, 6] < 300), :]
    ch16 = low_res[(low_res[:, 6] >= 300) & (low_res[:, 6] < 310), :]
    ch17 = low_res[(low_res[:, 6] >= 310) & (low_res[:, 6] < 320), :]
    ch18 = low_res[(low_res[:, 6] >= 320) & (low_res[:, 6] < 330), :]
    ch19 = low_res[(low_res[:, 6] >= 330) & (low_res[:, 6] < 340), :]
    ch20 = low_res[(low_res[:, 6] >= 340) & (low_res[:, 6] < 350), :]

    bits = np.unpackbits(truc[:, :57])
    show(bits.reshape(-1, 57 * 8)[::149, :], stripes=12)

    print "waiting..."
    raw_input()

    from trollimage.image import Image

    img = Image((ch14[:, 45:], ch12[:, 38:-7], ch09[:, 45:]), mode="RGB")
    img.stretch("logarithmic")
    img.show()
	#!/usr/bin/env python
	# -- coding: utf-8 --

	# Copyright (c) 2015 Martin Raspaud

	# Author(s):

	# Martin Raspaud <martin.raspaud@smhi.se>

	# This program is free software: you can redistribute it and/or modify
	# it under the terms of the GNU General Public License as published by
	# the Free Software Foundation, either version 3 of the License, or
	# (at your option) any later version.

	# This program is distributed in the hope that it will be useful,
	# but WITHOUT ANY WARRANTY; without even the implied warranty of
	# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	# GNU General Public License for more details.

	# You should have received a copy of the GNU General Public License
	# along with this program. If not, see <http://www.gnu.org/licenses/>.

	"""Read MERSI data from framed raw data.
	"""


	from datetime import datetime, timedelta
	import numpy as np
	import sys
	from reed import rs
	from trollimage.image import Image
	from bitstring import BitArray, BitStream
	import logging

	logger = None

	#from zfec import Decoder
	#from trollcast.formats.reedsolomon import *
	#from reedsolo import RSCodec

	table = np.array([0xff, 0x48, 0x0e, 0xc0, 0x9a, 0x0d, 0x70, 0xbc, 0x8e, 0x2c,
	0x93, 0xad, 0xa7, 0xb7, 0x46, 0xce, 0x5a, 0x97, 0x7d, 0xcc,
	0x32, 0xa2, 0xbf, 0x3e, 0x0a, 0x10, 0xf1, 0x88, 0x94, 0xcd,
	0xea, 0xb1, 0xfe, 0x90, 0x1d, 0x81, 0x34, 0x1a, 0xe1, 0x79,
	0x1c, 0x59, 0x27, 0x5b, 0x4f, 0x6e, 0x8d, 0x9c, 0xb5, 0x2e,
	0xfb, 0x98, 0x65, 0x45, 0x7e, 0x7c, 0x14, 0x21, 0xe3, 0x11,
	0x29, 0x9b, 0xd5, 0x63, 0xfd, 0x20, 0x3b, 0x02, 0x68, 0x35,
	0xc2, 0xf2, 0x38, 0xb2, 0x4e, 0xb6, 0x9e, 0xdd, 0x1b, 0x39,
	0x6a, 0x5d, 0xf7, 0x30, 0xca, 0x8a, 0xfc, 0xf8, 0x28, 0x43,
	0xc6, 0x22, 0x53, 0x37, 0xaa, 0xc7, 0xfa, 0x40, 0x76, 0x04,
	0xd0, 0x6b, 0x85, 0xe4, 0x71, 0x64, 0x9d, 0x6d, 0x3d, 0xba,
	0x36, 0x72, 0xd4, 0xbb, 0xee, 0x61, 0x95, 0x15, 0xf9, 0xf0,
	0x50, 0x87, 0x8c, 0x44, 0xa6, 0x6f, 0x55, 0x8f, 0xf4, 0x80,
	0xec, 0x09, 0xa0, 0xd7, 0x0b, 0xc8, 0xe2, 0xc9, 0x3a, 0xda,
	0x7b, 0x74, 0x6c, 0xe5, 0xa9, 0x77, 0xdc, 0xc3, 0x2a, 0x2b,
	0xf3, 0xe0, 0xa1, 0x0f, 0x18, 0x89, 0x4c, 0xde, 0xab, 0x1f,
	0xe9, 0x01, 0xd8, 0x13, 0x41, 0xae, 0x17, 0x91, 0xc5, 0x92,
	0x75, 0xb4, 0xf6, 0xe8, 0xd9, 0xcb, 0x52, 0xef, 0xb9, 0x86,
	0x54, 0x57, 0xe7, 0xc1, 0x42, 0x1e, 0x31, 0x12, 0x99, 0xbd,
	0x56, 0x3f, 0xd2, 0x03, 0xb0, 0x26, 0x83, 0x5c, 0x2f, 0x23,
	0x8b, 0x24, 0xeb, 0x69, 0xed, 0xd1, 0xb3, 0x96, 0xa5, 0xdf,
	0x73, 0x0c, 0xa8, 0xaf, 0xcf, 0x82, 0x84, 0x3c, 0x62, 0x25,
	0x33, 0x7a, 0xac, 0x7f, 0xa4, 0x07, 0x60, 0x4d, 0x06, 0xb8,
	0x5e, 0x47, 0x16, 0x49, 0xd6, 0xd3, 0xdb, 0xa3, 0x67, 0x2d,
	0x4b, 0xbe, 0xe6, 0x19, 0x51, 0x5f, 0x9f, 0x05, 0x08, 0x78,
	0xc4, 0x4a, 0x66, 0xf5, 0x58], dtype=np.uint8)


	def decode_str(arr):
	decoder = np.tile(table, np.ceil(len(arr) / 255.0))
	return (np.fromstring(arr, np.uint8) ^ decoder[:len(arr)]).tostring()


	def show(data, filename=None, stripes=0):
	"""Show the stetched data.
	"""

	norm_arr = np.array((data - data.min()) * 255.0 /
	(data.max() - data.min())).astype(np.uint8)
	if stripes:
	r_norm_arr = np.array(norm_arr)
	r_norm_arr[:, ::stripes] = np.uint8(255 - r_norm_arr[:, ::stripes])
	img = Image((r_norm_arr, norm_arr, norm_arr), mode="RGB")
	img.stretch("crude")
	else:
	from PIL import Image as pil
	img = pil.fromarray(norm_arr)

	#img = Image((ch14, ch12, ch09), mode="RGB")
	#img = Image((ch14[:, 38:], ch12[:, 38:], ch09[:, 38:]), mode="RGB")

	if filename:
	img.save(filename)
	else:
	img.show()


	def dec12to16(data):
	arr12 = data.astype(np.uint16).flat
	new_shape = list(data.shape[:-1]) + [(data.shape[-1] * 8) / 12]
	arr16 = np.zeros(new_shape, dtype=np.uint16)
	arr16.flat[::2] = np.left_shift(arr12[::3], 4) + \
	np.right_shift((arr12[1::3]), 4)
	arr16.flat[1::2] = np.left_shift((arr12[1::3] & 15), 8) + arr12[2::3]
	return arr16


	#############################################################

	from datetime import datetime, timedelta
	import numpy as np
	import sys
	#from reed import rs
	#from zfec import Decoder
	#from trollcast.formats.reedsolomon import *
	#from reedsolo import RSCodec, ReedSolomonError


	vcdu_type = np.dtype([('version', '>u2'),
	('vcdu count', '>u1', (3, )),
	('replay', '>u1'),
	('insert', '>u2'),
	('mpdu_hdr', '>u2')])

	ccsds_fy_hdr = np.dtype([('ccsds_version', '>u2'),
	('sequence', '>u2'),
	('packet_length', '>u2'),
	('days', '>u2'),
	('milliseconds', '>u4')])


	def cadu_info(line):
	print "=" * 3, "Cadu Info", "=" * 3
	print line["vcdu count"][0] * 256 * 256 + line["vcdu count"][1] * 256 + line["vcdu count"][2]
	print "vcid", line["version"] & (2 ** 6 - 1)
	print "spacecraft_id", (line["version"] >> 6) & (2 ** 8 - 1)
	print "version", (line["version"] >> 14) & 3
	if line["insert"]:
	print "Data crypted"
	else:
	print "Data open"
	print "offset", line["mpdu_hdr"] & (2 ** 11 - 1)


	def read_vcdu_hdr(packets):
	for packet in packets:
	line = np.fromstring(packet, vcdu_type, count=1)[0], packet[-882:]
	# cadu_info(line[0])
	# raw_input()
	yield line


	def reed_decode(packets):
	codec = rs.RSCoder(255, 223)
	res = bytearray("\0" * 892)
	#codec = RSCodec(32)

	for packet in packets:
	# do the decoding!
	yield packet[:892]
	# for i in range(4):
	# #print i
	# try:
	# print len(packet[i::4])
	# print codec.verify(packet[i::4])
	# decoded = codec.decode(packet[i::4])
	# if not decoded:
	# raise ValueError
	# #print decoded
	# res[i::4] = decoded
	# #print "decoded"
	# #raw_input()
	# except IOError:
	# print "uncorrectable"
	# res[i::4] = packet[i:892:4]
	# except ReedSolomonError as err:
	# #print "uncorrectable:", err
	# res[i::4] = packet[i:892:4]
	# print str(res[:892]) == packet[:892]
	# yield str(res[:892])


	def pn_decode(packets, dlen=0):
	for packet in packets:
	if dlen < 0:
	yield packet
	else:
	yield packet[:dlen] + decode_str(packet[dlen:])


	def frame_sync(packets):
	"""Phony frame_sync.
	"""
	for packet in packets:
	yield packet[4:]


	def file_reader(filename, plen=1024):
	with open(filename) as fd:
	while True:
	packet = fd.read(plen)
	if packet:
	yield packet
	else:
	return


	def sort_vc(packets):
	channels = {}
	cnt = 0
	for hdr, packet in packets:
	vcid = hdr["version"] & (2 ** 6 - 1)
	channels.setdefault(vcid, []).append((hdr, packet))
	if cnt % 1000 == 0:
	logger.debug(str(cnt))
	cnt += 1
	return channels


	def gen_strips(data):
	frame_sync = "\xff\xe0\x01\xff\xe0\x01\xff\xe0\x01"

	bla = BitStream(bytes=data[:55000])
	offset = bla.find(BitArray(bytes=frame_sync))[0] % 8
	logger.debug("Bit offset %d", offset)
	if offset != 0:
	array = np.fromstring(data, dtype=np.uint8)
	data = (((array[:-1] & (2**(8 - offset) - 1)) << offset) + (array[1:] >> (8 - offset))).astype(np.uint8).tostring()
	prev = 0
	while True:
	try:
	idx = data.index(frame_sync, prev + len(frame_sync))
	except ValueError:
	break
	yield data[prev:idx]
	prev = idx


	def read5(filename, plen, pn_start):
	channels = sort_vc(
	read_vcdu_hdr(reed_decode(pn_decode(frame_sync(file_reader(filename, plen=plen)), pn_start))))
	logger.debug('Channel keys: %s', str(channels.keys()))

	return channels

	if __name__ == '__main__':

	import argparse

	parser = argparse.ArgumentParser()
	parser.add_argument("-c", "--cadu-len", help="Length of the cadu package in bytes (usually 1024 bytes)",
	type=int, default=1024)
	parser.add_argument("-p", "--pseudo-noise-start", help="Start of the pseudo-noise chunk, usually 0. Default is -1 for no decoding.",
	type=int, default=-1)
	parser.add_argument("-v", "--verbose", help="print debug messages too",
	action="store_true")
	parser.add_argument("filename", help="name of the file to read")
	opts = parser.parse_args()

	handler = logging.StreamHandler()
	handler.setFormatter(logging.Formatter("[%(levelname)s: %(asctime)s :"
	" %(name)s] %(message)s",
	'%Y-%m-%d %H:%M:%S'))
	if opts.verbose:
	loglevel = logging.DEBUG
	else:
	loglevel = logging.INFO
	handler.setLevel(loglevel)
	logging.getLogger('').setLevel(loglevel)
	logging.getLogger('').addHandler(handler)
	logger = logging.getLogger("fy3_mersi_reader")


	chs = read5(opts.filename, opts.cadu_len, opts.pseudo_noise_start)

	logger.debug(str(chs.keys()))

	data = "".join(dat[1] for dat in chs[3])

	res = []
	low_res = []
	truc = []
	rest = []
	truc_str = []
	logger.info('Stripping data')
	for strip in gen_strips(data):
	if len(strip) == 12345: #12345:
	a = np.fromstring(strip, dtype=np.uint8)
	res.append(dec12to16(a))

	elif len(strip) == 3129:
	a = np.fromstring(strip, dtype=np.uint8)
	low_res.append(dec12to16(a))
	truc.append(a)
	elif len(strip) == 54531:
	a = np.fromstring(strip[:3129], dtype=np.uint8)
	low_res.append(dec12to16(a))
	truc_str.append(strip[3129:])
	rest.append(dec12to16(np.fromstring(strip[3129:], dtype=np.uint8)))
	truc = np.vstack(truc)
	low_res = np.vstack(low_res)
	res = np.vstack(res)
	restarr = np.vstack(rest)

	#for i in range(100):
	#a = truc[150 * i, 12:30].astype(np.uint32)
	#a = np.take(a, np.arange(6) + 3, mode="wrap")
	# milli = ((a[3] & 15) << 28) + ((a[4] >> 4) << 24) + \
	# ((a[1] & 15) << 20) + (a[2] << 12) + (a[0] << 4) + a[1] >> 4
	#milli = (((((a[3] >> 4) << 12) + a[3]) << 12) + a[2])

	#milli = 0
	#day = ((a[4] & 15) << 12) + (a[5] << 4) + (a[3] >> 4)

	# print a
	# half = (((a[:-1] & 15) << 4) + (a[1:] >> 4)).astype(np.uint8)
	# print half
	# print dec12to16(a)
	# day = (half[4] << 8) + half[5]
	# milli = 0

	#day = (a[0] << 8) + a[1]
	#milli = (a[2] << 24) + (a[3] << 16) + (a[4] << 8) + a[5]

	#thedate = datetime(
	# 2000, 1, 1) + timedelta(days=int(day) - 1, milliseconds=int(milli))
	#print thedate

	#print "waiting..."
	#raw_input()

	res_chn = ["01", "02", "03", "04", "05"]
	low_res_chn = ["01", "02", "03", "04", "05"]

	ch01 = res[res[:, 6] < 40, :]
	ch02 = res[(res[:, 6] >= 40) & (res[:, 6] < 80), :]
	ch03 = res[(res[:, 6] >= 80) & (res[:, 6] < 120), :]
	ch04 = res[(res[:, 6] >= 120) & (res[:, 6] < 160), :]
	ch05 = res[(res[:, 6] >= 160) & (res[:, 6] < 200), :]

	ch06 = low_res[low_res[:, 6] < 210, :]
	ch07 = low_res[(low_res[:, 6] >= 210) & (low_res[:, 6] < 220), :]
	ch08 = low_res[(low_res[:, 6] >= 220) & (low_res[:, 6] < 230), :]
	ch09 = low_res[(low_res[:, 6] >= 230) & (low_res[:, 6] < 240), :]
	ch10 = low_res[(low_res[:, 6] >= 240) & (low_res[:, 6] < 250), :]
	ch11 = low_res[(low_res[:, 6] >= 250) & (low_res[:, 6] < 260), :]
	ch12 = low_res[(low_res[:, 6] >= 260) & (low_res[:, 6] < 270), :]
	ch13 = low_res[(low_res[:, 6] >= 270) & (low_res[:, 6] < 280), :]
	ch14 = low_res[(low_res[:, 6] >= 280) & (low_res[:, 6] < 290), :]
	ch15 = low_res[(low_res[:, 6] >= 290) & (low_res[:, 6] < 300), :]
	ch16 = low_res[(low_res[:, 6] >= 300) & (low_res[:, 6] < 310), :]
	ch17 = low_res[(low_res[:, 6] >= 310) & (low_res[:, 6] < 320), :]
	ch18 = low_res[(low_res[:, 6] >= 320) & (low_res[:, 6] < 330), :]
	ch19 = low_res[(low_res[:, 6] >= 330) & (low_res[:, 6] < 340), :]
	ch20 = low_res[(low_res[:, 6] >= 340) & (low_res[:, 6] < 350), :]

	bits = np.unpackbits(truc[:, :57])
	show(bits.reshape(-1, 57 * 8)[::149, :], stripes=12)

	print "waiting..."
	raw_input()

	from trollimage.image import Image

	img = Image((ch14[:, 45:], ch12[:, 38:-7], ch09[:, 45:]), mode="RGB")
	img.stretch("logarithmic")
	img.show()