misc.ecl

   1"""Read Schlumberger Eclipse output files.
   2"""
   3from __future__ import division
   4import argparse
   5import codecs
   6import collections
   7import datetime
   8import fnmatch
   9import logging
  10import numpy
  11import numpy.ma
  12import os
  13import os.path as path
  14import struct
  15import sys
  16
  17
  18# import these from the common package, since we may need to use the same
  19# enumeration across file formats
  20from .ecl_common import Phase, Prop
  21
  22
  23# add a valid log in case we are not run through the main program which
  24# sets up one for us
  25log = logging.getLogger(__name__)  # pylint: disable=invalid-name
  26log.addHandler(logging.NullHandler())
  27
  28
  29# each keyword record in the file is associated with a name and an index
  30# (since a keyword may occur more than once)
  31_DataKey = collections.namedtuple('_DataKey', [
  32    'kwd',  # keyword present at this position
  33    'seq',  # sequence number (occurrence in file)
  34])
  35
  36
  37# for every (indexed) keyword, we want to know where to read the data
  38# in the file, as well as how much data there is
  39_DataDescr = collections.namedtuple('_DataDescr', [
  40    'pos',  # file position of this record
  41    'num',  # number of elements
  42    'typ',  # type of data
  43    'rcs',  # number of records containing data
  44])
  45
  46
  47def _read_descr(fileobj):
  48    """Read the the opened file and build a descriptor record of the
  49    following data array.
  50
  51    :param fileobj: File object opened in binary reading mode.
  52    :returns:       Tuple of keyword, file position, number of items,
  53                    and their type.
  54    :rtype:         (str, int, int, str)
  55    """
  56    # read the length of this record, as a signed quad; if we encounter
  57    # end of file, then return None as the sentinel descriptor
  58    quad = bytes(fileobj.read(4))
  59    if len(quad) == 0:
  60        return None, None, None, None
  61    rec_len, = struct.unpack('>i', quad)
  62
  63    # sub-records should not occur
  64    assert (rec_len == 16)
  65
  66    # read the next 16 bytes as descriptor
  67    kwd, = struct.unpack('8s', fileobj.read(8))
  68    num, = struct.unpack('>i', fileobj.read(4))
  69    typ, = struct.unpack('4s', fileobj.read(4))
  70
  71    # for Python 3 compatibility
  72    kwd = codecs.ascii_decode(kwd)[0]
  73    typ = codecs.ascii_decode(typ)[0]
  74
  75    # seek forward to the end of the record
  76    fileobj.seek(abs(rec_len) - 16, 1)
  77
  78    # read the trailing length as well
  79    trail, = struct.unpack('>i', fileobj.read(4))
  80
  81    # verify that the trailing length is the same as the
  82    assert (trail == rec_len)
  83
  84    # get the position of the next record in the file; this is where
  85    # we want to seek to, to start reading the data array
  86    pos = fileobj.tell()
  87
  88    # return the descriptor of this array
  89    return (kwd, pos, num, typ)
  90
  91
  92# for each type of data, we need to know how much to read
  93_DataType = collections.namedtuple('_DataType', [
  94    'siz',  # number of bytes
  95    'fmt',  # struct formatting string
  96    'nch',  # not character type, which needs special processing
  97    'dsk',  # NumPy data type when stored on disk (big-endian)
  98    'mem',  # corresponding NumPy data type in memory
  99])
 100
 101
 102# Python 3 uses Unicode for characters
 103_STR_TYPE = 'S' if sys.hexversion < 0x03000000 else 'U'
 104
 105
 106# actual data types found in Eclipse files
 107_TYPES = {
 108    'INTE': _DataType(siz=4, fmt='i', nch=True,
 109                      dsk=numpy.dtype('>i4'), mem=numpy.int32),
 110    'REAL': _DataType(siz=4, fmt='f', nch=True,
 111                      dsk=numpy.dtype('>f4'), mem=numpy.float32),
 112    'LOGI': _DataType(siz=4, fmt='i', nch=True,
 113                      dsk=numpy.dtype('>i4'), mem=numpy.int32),
 114    'DOUB': _DataType(siz=8, fmt='d', nch=True,
 115                      dsk=numpy.dtype('>f8'), mem=numpy.float64),
 116    'CHAR': _DataType(siz=8, fmt='c', nch=False,
 117                      dsk=numpy.dtype('|S8'), mem=(
 118                          numpy.dtype('<' + _STR_TYPE + '8'))),
 119    'MESS': _DataType(siz=0, fmt='x', nch=True,
 120                      dsk=numpy.void, mem=numpy.void),
 121}
 122
 123
 124def _skip_rec(fileobj, num, typ):
 125    """Skip one or more records for a certain number of entries.
 126
 127    :param fileobj: File object opened in binary read mode.
 128    :param num:     Number of items that should be skipped.
 129    :param typ:     Data type of the entries.
 130    :returns:       Number of data records that were skipped.
 131    """
 132    # total number of bytes to be skipped
 133    remaining = num * _TYPES[typ].siz
 134
 135    # number of records these bytes were stored in
 136    rcs = 0
 137
 138    # read until we have processed the entire data array
 139    while remaining > 0:
 140        # read the number of bytes in this record
 141        rec_len, = struct.unpack('>i', fileobj.read(4))
 142
 143        # this part of the data array is now read
 144        remaining -= rec_len
 145        rcs += 1
 146
 147        # actually skip those bytes, plus the trailing record length
 148        fileobj.seek(rec_len + 4, os.SEEK_CUR)
 149
 150    # we should not end up in the middle of a record!
 151    assert (remaining == 0)
 152
 153    # store the number of records so that we know exactly how many
 154    # records to read next time
 155    return rcs
 156
 157
 158def _build_cat(fileobj):
 159    """Build a catalog over where in the data file each record is.
 160
 161    :param fileobj:  File object opened in binary read mode.
 162    :returns:        Tuple of dictionaries with the number of occurrences
 163    :rtype:          (dict [str] -> int, dict [_DataKey] -> _DataDescr)
 164    """
 165    cnt = {}  # number of types we have seen each keyword
 166    cat = {}  # where is the file each keyword is located
 167    while True:
 168        # read the record containing the data descriptor first
 169        kwd, pos, num, typ = _read_descr(fileobj)
 170
 171        # stop once we have reached the end of the file (no more records)
 172        if kwd is None:
 173            break
 174
 175        # skip over the data records for this key; we are only interested
 176        # in building the index at this stage
 177        rcs = _skip_rec(fileobj, num, typ)
 178
 179        # update the number of occurrences of this keyword
 180        cnt[kwd] = cnt[kwd] + 1 if kwd in cnt else 1
 181
 182        # create a descriptor for this record
 183        cat[_DataKey(kwd=kwd, seq=cnt[kwd] - 1)] = (
 184            _DataDescr(pos=pos, num=num, typ=typ, rcs=rcs))
 185
 186    return (cnt, cat)
 187
 188
 189def _read_rec(fileobj, descr):
 190    """Read a set of records for a descriptor.
 191
 192    :param fileobj:  File object opened in binary read mode.
 193    :param descr:    Descriptor for where the property is stored.
 194    :type descr:     _DataDescr
 195    """
 196    # seek to the appropriate place in the file
 197    fileobj.seek(descr.pos, os.SEEK_SET)
 198
 199    # get a description of the type of data
 200    rec_typ = _TYPES[descr.typ]
 201
 202    # pre-allocate an array so we can just read straight into it
 203    data = numpy.empty(descr.num, dtype=rec_typ.dsk)
 204
 205    # start reading into the front of the array
 206    fst = 0
 207
 208    # read as many records as needed from here on
 209    for rec_ndx in range(descr.rcs):  # pylint: disable=unused-variable
 210        # how many bytes are there in this particular record
 211        rec_siz, = struct.unpack('>i', fileobj.read(4))
 212
 213        # an item cannot be split over two physical records
 214        assert (rec_siz % rec_typ.siz == 0)
 215
 216        # how many *items* are there in this record
 217        rec_num = int(abs(rec_siz) / rec_typ.siz)
 218
 219        # now we know which range we are going to read into;
 220        # bulk load it directly from file. character records need special
 221        # handling.
 222        if rec_typ.nch:
 223            data[fst:(fst + rec_num)] = (
 224                numpy.fromfile(fileobj, dtype=rec_typ.dsk, count=rec_num))
 225        else:
 226            # load the characters as uint8 raw bytes from the file
 227            raw = numpy.fromfile(fileobj, dtype=numpy.uint8,
 228                                 count=rec_num * rec_typ.siz)
 229            for item_ndx in range(rec_num):
 230                # start and end iterator for where in the raw array the bytes
 231                # for this particular item in the array is; each of the strings
 232                # are considered to be *exactly* 8 bytes by Eclipse
 233                beg = (item_ndx + 0) * rec_typ.siz
 234                end = (item_ndx + 1) * rec_typ.siz
 235
 236                # convert these raw characters into a sensible string, and
 237                # store them in the output array
 238                str_val = codecs.ascii_decode(bytearray(raw[beg: end]))[0]
 239                data[fst + item_ndx] = str_val.rstrip()
 240
 241        # skip the trailing record size
 242        fileobj.seek(4, os.SEEK_CUR)
 243
 244        # move to the next window
 245        fst += rec_num
 246
 247    # if we loaded big-endian data on a little-endian machine,
 248    # swap the entire array in-place now, and then cast it
 249    # (with a view) to the native type so we see the bits as
 250    # the right values again
 251    if rec_typ.nch:
 252        if sys.byteorder != 'big':
 253            data.byteswap(True)
 254        data = data.view(rec_typ.mem)
 255    else:
 256        # character data cannot be casted, but must be converted into the
 257        # correct number of bytes per character
 258        data = data.astype(rec_typ.mem)
 259
 260    return data
 261
 262
 263class EclipseFile (object):
 264    """Low-level class to read records from binary files.
 265
 266    Access to this object must be within a monitor (`with`-statement).
 267    """
 268
 269    def __init__(self, root, ext):
 270        """Initialize file object from a path in the filesystem.
 271
 272        :param root:  Stem of the file name (including directory)
 273        :param ext:   Extension of the file to read from.
 274        """
 275        # keep the file open (and locked) to read from it
 276        self.filename = '{0}.{1}'.format(root, ext.upper())
 277        self.fileobj = open(self.filename, 'rb')
 278
 279        # index the file so that we can find properties easily
 280        log.debug("Indexing data file \"%s\"", self.filename)
 281        self.cnt, self.cat = _build_cat(self.fileobj)
 282
 283    def __enter__(self):
 284        self.fileobj.__enter__()
 285        return self
 286
 287    def __exit__(self, typ, val, traceback):
 288        """Close the underlaying file object when we go out of scope."""
 289        self.fileobj.__exit__(typ, val, traceback)
 290
 291    def get(self, kwd, seq=0):
 292        """Read a (potentially indexed) keyword from file.
 293
 294        :param kwd:  Keyword.
 295        :type kwd:   str
 296        :param seq:  Sequence number, starting from zero.
 297        :type seq:   int
 298
 299        >>> with EclipseFile ('foo', 'EGRID') as grid:
 300        >>>     zcorn = grid.get ('ZCORN')
 301        """
 302        # make sure that the keyword is exactly 8 chars
 303        kwd = "{0:<8s}".format(kwd[:8].upper())
 304
 305        # locate the data in the file and read them there, returning
 306        # None if we didn't find any matching record in the catalog
 307        key = _DataKey(kwd=kwd, seq=seq)
 308        if key in self.cat:
 309            descr = self.cat[key]
 310            data = _read_rec(self.fileobj, descr)
 311            return data
 312        else:
 313            return None
 314
 315    def dump(self, positional=False, fileobj=sys.stdout):
 316        """Dump catalog contents of records in the datafile.
 317
 318        :param positional:  True if the keywords should be sorted
 319                            on position in the file.
 320        :type positional:   bool
 321
 322        >>> f = EclipseFile ("foo", 'INIT')
 323        >>> f.dump ()
 324        """
 325        if positional:
 326            # build a dictionary of positions and they record that
 327            # are at each of those positions. by using the position
 328            # as the key of the dictionary, we can sort them and then
 329            # retrive the keys in positionally sorted order
 330            by_pos = {}
 331            for key in self.cat.keys():
 332                descr = self.cat[key]
 333                by_pos[descr.pos] = key
 334            key_list = [by_pos[pos] for pos in sorted(by_pos.keys())]
 335        else:
 336            # sort the list of keywords alphabetically and then
 337            # by sequence number
 338            key_list = sorted(self.cat.keys())
 339
 340        # write each entry in the catalog to the console
 341        for key in key_list:
 342            kwd, seq = key
 343            descr = self.cat[key]
 344            total = self.cnt[kwd]
 345
 346            # print different text if there is only one of this keyword
 347            if total == 1:
 348                fmt = '{0:<8s}  {5:7s}  {2:>8d} * {3:4s}\n'
 349            else:
 350                fmt = '{0:<8s}  {1:03d}/{4:03d}  {2:>8d} * {3:4s}\n'
 351
 352            # write information for this record
 353            fileobj.write(fmt.format(
 354                kwd, seq + 1, descr.num, descr.typ, total, ''))
 355
 356
 357# pylint: disable=too-few-public-methods, too-many-instance-attributes
 358class EclipseGrid (object):
 359    """Corner-point geometry data from an Eclipse Extensive Grid file."""
 360
 361    def __init__(self, root):
 362        # save the root variable for later
 363        self.root = root
 364
 365        # read the extent of the grid right away since we need that for
 366        # later (when reading properties).
 367        with EclipseFile(root, "EGRID") as egrid:
 368            # read the grid header
 369            grid_head = egrid.get('GRIDHEAD')
 370
 371            # only corner-point grids are supported
 372            assert (grid_head[0] == 1)
 373
 374            # get dimensions of the grid
 375            self.ni = grid_head[1]  # pylint: disable=invalid-name
 376            self.nj = grid_head[2]  # pylint: disable=invalid-name
 377            self.nk = grid_head[3]  # pylint: disable=invalid-name
 378            log.info("Grid dimension is %d x %d x %d",
 379                     self.ni, self.nj, self.nk)
 380
 381            # also store a shape tuple which describes the grid cube
 382            self.shape = (self.nk, self.nj, self.ni)
 383
 384            # read the structure of the grid
 385            log.debug("Reading pillar coordinates")
 386            coord = egrid.get('COORD')
 387            self.coord = numpy.reshape(coord,
 388                                       (self.nj + 1, self.ni + 1, 2, 3))
 389
 390            log.debug("Reading corner depths")
 391            zcorn = egrid.get('ZCORN')
 392            self.zcorn = numpy.reshape(zcorn,
 393                                       (self.nk, 2, self.nj, 2, self.ni, 2))
 394
 395            log.debug("Reading active flag")
 396            actnum = egrid.get('ACTNUM')
 397
 398            # make it a true flag, so we can use it for binary indexing
 399            actnum = actnum.astype(bool)
 400
 401            # now we keep the actnum array in its proper shape,
 402            # corresponding to the grid cube
 403            self.actnum = numpy.reshape(actnum, self.shape)
 404
 405            # create a common mask in the grid that we can use in
 406            # masked arrays for other properties.
 407            self.mask = numpy.logical_not(self.actnum)
 408
 409            # restart properties are only saved for the active elements,
 410            # so we can cache this number to compare
 411            self.num_active = numpy.sum(self.actnum)
 412            log.info("Grid has %d active cells", self.num_active)
 413
 414    def grid(self):
 415        """\
 416        Create a grid structure from the information read from file.
 417
 418        :return:  Grid structure
 419        :rtype:   dict
 420
 421        >>> # convert from .EGRID to .grdecl:
 422        >>> import pyresito.io.ecl as ecl
 423        >>> import pyresito.io.grdecl as grdecl
 424        >>> grdecl.write('foo.grdecl', ecl.EclipseGrid('FOO').grid())
 425        """
 426        # we have all the information that is required in the correct
 427        # format already as soon as we read it from file
 428        dimens = numpy.array([self.shape[2],
 429                              self.shape[1],
 430                              self.shape[0]], dtype=numpy.int32)
 431
 432        # return grid in the same format as if it was read directly from
 433        # the original corner-point grid definition (.grdecl)
 434        return {'DIMENS': dimens,
 435                'COORD':  self.coord.copy(),
 436                'ZCORN':  self.zcorn.copy(),
 437                'ACTNUM': self.actnum.copy(),
 438                }
 439
 440
 441# Eclipse 300 on Windows initialize the intehead array to this value
 442# and doesn't bother to write all values
 443_UNINITIALIZED = -2345
 444
 445
 446def _intehead_date(intehead):
 447    """Get the date coded in the file header."""
 448    year = intehead[66]
 449    month = intehead[65]
 450    day = intehead[64]
 451
 452    # libecl used in OPM seems to write only date information, not time,
 453    # unless it is necessary; the header will then be only 95 integers long
 454    # see function ecl_init_file_alloc_INTEHEAD at lines 56-59 in file
 455    # lib/ecl/ecl_init_file.cpp in libecl release 2018.10 commit 77ca7d5
 456    if (len(intehead) > 208):
 457        hour = intehead[206]
 458        minute = intehead[207]
 459
 460        # take into account special values that indicate uninitialized field
 461        hour = 0 if hour == _UNINITIALIZED else hour
 462        minute = 0 if minute == _UNINITIALIZED else minute
 463    else:
 464        hour = 0
 465        minute = 0
 466
 467    # init files from the Windows version doesn't have this field (?)
 468    if (len(intehead) >= 411):
 469        msecs = intehead[410]
 470    else:
 471        msecs = 0
 472
 473    # split microseconds into seconds and remaining microseconds
 474    secs = int(msecs / 1e6)
 475
 476    # spurious values found in this field which makes it unaligned with
 477    # the start date; better just drop the microsecond field
 478    msecs = 0   # msecs -= int (secs * 1e6)
 479
 480    return datetime.datetime(year, month, day,
 481                             hour, minute, secs, msecs)
 482
 483
 484# pylint: disable=too-few-public-methods, too-many-instance-attributes
 485class EclipseData (object):
 486    """Base class for both static and recurrent data."""
 487
 488    def __init__(self, grid, root, ext):
 489        self.grid = grid  # grid extent information
 490        self.root = root  # underlaying file name
 491        self.ext = ext    # connected file with extent of grid
 492
 493        # list of components isn't loaded yet (and may never be)
 494        self.comp = None
 495
 496    def components(self):
 497        """Components that exist in the restart file. Components used in the
 498        simulation are stored in all the restart files instead of once in the
 499        init file, since it is the restart file that hold component fields.
 500        """
 501        # delay-load the list of components; if it already is loaded,
 502        # then just return the stored one. we cache component names since
 503        # they are needed to lookup property names.
 504        if self.comp is None:
 505            # load that datafile
 506            with EclipseFile(self.root, self.ext) as store:
 507                raw_comp = store.get('ZCOMPS')
 508
 509            # create a dictionary over all the components, in the order they
 510            # were specified in the file, since component-based properties have
 511            # this index as part of their name.
 512            self.comp = collections.OrderedDict(
 513                [(str(nom), ndx + 1) for ndx, nom in enumerate(raw_comp)])
 514
 515        return self.comp
 516
 517    def _comp_phase(self, selector):  # pylint: disable=no-self-use
 518        """Code indicating the phase, when we are asking for a component-based
 519        property. This code is different than when asking for a phase-based
 520        property.
 521        """
 522        # second argument is the component, (optional) third argument is
 523        # the phase. when it comes to listing mole fractions, simulators (for
 524        # some unknown reason) uses other mnemonics than O, W, G to designate
 525        # phases.
 526        if len(selector) == 2:
 527            mph = 'Z'
 528        elif selector[2] == Phase.oil:
 529            mph = 'X'
 530        elif selector[2] == Phase.wat:
 531            mph = 'A'
 532        elif selector[2] == Phase.gas:
 533            mph = 'Y'
 534        else:
 535            assert False
 536
 537        return mph
 538
 539    def _get_prop_name(self, selector):  # pylint: disable=no-self-use
 540        """Compose an internal name for the property as specified by a
 541        high-level selector.
 542
 543        :param selector:  Selector tuple, e.g. (Prop.mole, 'CO2', Phase.gas)
 544        :param names:     Dictionary of defined names in the case. There must
 545                          be an entry "components" containing the names of the
 546                          components in the case files.
 547        :returns:         Name of the property to be queried, e.g. "ZMF2"
 548        """
 549        # if the selector is a tuple, then use the first member to get a
 550        # function that can parse the others
 551        if isinstance(selector, tuple):
 552            # take the full selector (and the name dictionary to supplement)
 553            # and compose the final name. the selector can therefore be thought
 554            # of as an s-expr which is evaluated here
 555            if selector[0] == Prop.pres:
 556                name = 'PRESSURE'
 557            elif selector[0] == Prop.sat:
 558                # second argument is the phase
 559                name = 'S' + selector[1]
 560            elif selector[1] == Prop.mole:
 561                # get the phase scope of the component property
 562                scope = self._comp_phase(selector)
 563
 564                # use dictionary to find the component index number
 565                ndx = ('' if selector[1] is None
 566                       else str(self.components()[selector[1]]))
 567
 568                name = scope + 'MF' + ndx
 569            else:
 570                assert False
 571
 572        else:
 573            # if it is not a selector, then it is (probably) a string, and
 574            # we forward it without any formatting, for compatibility
 575            name = selector
 576
 577        return name
 578
 579    def cell_data(self, selector):
 580        """Get a field property for every cell at this restart step.
 581
 582        :param selector: Specification of the property to be loaded. This is a
 583                         tuple starting with a Prop, and then some context-
 584                         dependent items.
 585        :returns:        Array of the data with inactive cells masked off.
 586
 587        >>> pres = case.cell_data ((Prop.pres,))
 588        >>> swat = case.cell_data ((Prop.sat, Phase.wat))
 589        >>> xco2 = case.cell_data ((Prop.mole, 'CO2', Phase.gas))
 590        """
 591        # get the internal name of the property requested
 592        propname = self._get_prop_name(selector)
 593
 594        # load the data itself
 595        with EclipseFile(self.root, self.ext) as store:
 596            active_data = store.get(propname)
 597
 598        # if there is a fully specified array, then just use the data
 599        if active_data.shape[0] == numpy.prod(self.grid.shape):
 600            data = numpy.reshape(active_data, self.grid.shape)
 601
 602        # if the data is "compressed", i.e. only the active elements are
 603        # stored, then uncompress it before storing it in a masked array
 604        elif active_data.shape[0] == self.grid.num_active:
 605            # allocate an array to hold the entire array. fill with zeros
 606            # since this value exists for all datatypes (nan doesn't)
 607            data = numpy.zeros(self.grid.shape, dtype=active_data.dtype)
 608
 609            # assign only the values that are actually active
 610            data[self.grid.actnum] = active_data
 611        else:
 612            assert (False)
 613
 614        # return an array where the non-active elements are masked away
 615        return numpy.ma.array(data=data, dtype=data.dtype,
 616                              mask=self.grid.mask)
 617
 618    def field_data(self, propname):
 619        """Get a property for the entire field at this restart step.
 620
 621        :param propname: Name of the property to be loaded.
 622        :returns:        Array of the data.
 623        """
 624        # load the data itself
 625        with EclipseFile(self.root, self.ext) as store:
 626            return store.get(propname)
 627
 628    def summary_data(self, propname):
 629        """Get a property from the summary file at this restart step.
 630
 631        :param propname: Name of the property to be loaded. This is on the form
 632                         'mnemonic well', e.g. 'WWIR I05'. Alternatively, one
 633                         propname can be either only well or only mnemonic.
 634                         Then the value for all mnememonics or all wells are
 635                         given, e.g., propname='WWIR' returns WWIR for all
 636                         wells.
 637        :returns:        Array of the data.
 638        """
 639        # need to find the mneumonics and well info from the specification file
 640        with EclipseFile(self.root, 'SMSPEC') as store:
 641            mnemonic = store.get('KEYWORDS')
 642            well = store.get('WGNAMES')
 643
 644        # split the propname
 645        prop_elem = propname.upper().split()
 646
 647        assert len(prop_elem) < 3
 648
 649        if len(prop_elem) == 1:  # only well or mnemonic is given
 650            if prop_elem in mnemonic:
 651                ind = numpy.where(mnemonic == prop_elem)
 652            elif prop_elem in well:
 653                ind = numpy.where(well == prop_elem)
 654
 655        elif len(prop_elem) == 2:
 656            # find all indexes for the mnemonic
 657            tmp_ind = numpy.where(mnemonic == prop_elem[0])
 658
 659            # fine all indexes for the well
 660            tmp_ind2 = numpy.where(well == prop_elem[1])
 661
 662            # only one element is equal
 663            ind = tmp_ind[0][numpy.array(
 664                [elem in tmp_ind2[0] for elem in tmp_ind[0]])]
 665
 666        # The challenge with the summary files is that it contains vectors of
 667        # values for all "ministeps" controlled by the internal time stepping
 668        # procedure. We are only interested in values at the report time, e.g.,
 669        # the final "ministep". It it, however, possible to collect all values,
 670        # e.g., for number of newton iterations.
 671        with EclipseFile(self.root, self.ext) as store:
 672            # Since we do not know the number of timesteps we must find how
 673            # many times the PARAM keyword has been given
 674            total = store.cnt['PARAMS  ']  # Gives number of ministeps
 675
 676            # Get the final instance, need to remove one since python counts
 677            # from zero
 678            vec_dat = store.get(kwd='PARAMS  ', seq=total-1)
 679
 680        return vec_dat[ind]
 681
 682
 683def _intehead_phases(intehead):
 684    """Get list of faces from the integer header"""
 685    # no phases discovered yet
 686    phs = []
 687
 688    # build list of phases from the flags in the header
 689    if intehead[14] & 1 == 1:  # bit 0 set if oil present
 690        phs.append(Phase.oil)
 691    if intehead[14] & 2 == 2:  # bit 1 set if water present
 692        phs.append(Phase.wat)
 693    if intehead[14] & 4 == 4:  # bit 2 set if gas present
 694        phs.append(Phase.gas)
 695
 696    return phs
 697
 698
 699# pylint: disable=too-few-public-methods, too-many-instance-attributes
 700class EclipseInit (EclipseData):
 701    """Read information from static data (init) file."""
 702
 703    def __init__(self, root):
 704        # bootstrap ourself; if the base class is reading a property
 705        # that require reformatting, then we have all the information
 706        # needed right here!
 707        super(EclipseInit, self).__init__(self, root, 'INIT')
 708
 709        # read properties from file (notice that file is closed afterwards)
 710        with EclipseFile(self.root, self.ext) as store:
 711            intehead = store.get('INTEHEAD')
 712            porv = store.get('PORV')
 713
 714        # get the dimensions from the grid
 715        self.ni = intehead[8]  # pylint: disable=invalid-name
 716        self.nj = intehead[9]  # pylint: disable=invalid-name
 717        self.nk = intehead[10]  # pylint: disable=invalid-name
 718
 719        # notice that we store the grid in C order, i.e. the dimensions
 720        # are swapped compared with what they are in Eclipse (Fortran)
 721        self.shape = (self.nk, self.nj, self.ni)
 722
 723        # get the starting date of the simulation
 724        self.start_date = _intehead_date(intehead)
 725
 726        # get the phases in the simulation
 727        self.phases = _intehead_phases(intehead)
 728
 729        # number of active elements is explicitly listed
 730        self.num_active = intehead[11]
 731
 732        # get the pore volume and reformat it into the shape of the grid
 733        self.pore_vol = numpy.reshape(porv, self.shape)
 734
 735        # create a mask and actnum array based on the pore volume
 736        # (active cells must have strictly positive volume)
 737        self.actnum = numpy.greater(self.pore_vol, 0.)
 738        self.mask = numpy.logical_not(self.actnum)
 739        assert (numpy.sum(self.actnum.ravel()) == self.num_active)
 740
 741
 742class EclipseRestart (EclipseData):
 743    """Read information from a recurrent data (restart) file."""
 744
 745    def __init__(self, grid, seq):
 746        """
 747        :param grid:  Initialization file which contains grid dimensions.
 748        :type grid:   EclipseInit (or EclipseGrid)
 749        :param seq:   Run number.
 750        :type seq:    int
 751        """
 752        super(EclipseRestart, self).__init__(grid, grid.root,
 753                                             "X{0:04d}".format(seq))
 754
 755    def date(self):
 756        """Simulation date the restart file is created for."""
 757        # dates are stored in the header
 758        with EclipseFile(self.root, self.ext) as store:
 759            intehead = store.get('INTEHEAD')
 760
 761        # convert Eclipse date field to a Python date object
 762        return _intehead_date(intehead)
 763
 764
 765class EclipseSummary (EclipseData):
 766    """Read information from a recurrent data (summary) file."""
 767
 768    def __init__(self, grid, seq):
 769        """
 770        :param grid:  Initialization file which contains grid dimensions.
 771        :type grid:   EclipseInit (or EclipseGrid)
 772        :param seq:   Run number.
 773        :type seq:    int
 774        """
 775        super(EclipseSummary, self).__init__(grid, grid.root,
 776                                             "S{0:04d}".format(seq))
 777
 778    def date(self):
 779        """Simulation date the restart file is created for."""
 780        # dates are stored in the header
 781        with EclipseFile(self.root, self.ext) as store:
 782            intehead = store.get('INTEHEAD')
 783
 784        # convert Eclipse date field to a Python date object
 785        return _intehead_date(intehead)
 786
 787
 788def _quick_date(fileobj):
 789    """Quick scan of a restart file to determine its date.
 790
 791    :param fileobj:  File object opened in binary read mode.
 792    """
 793    # INTEHEAD record is always placed first in the file; we have to read
 794    # four sectors; first sector contains the signature, second sector
 795    # (around 24 + 64*4 = 280) contains the day and (around 24 + 206*4 = 848)
 796    # time and fourth sector (around 24 + 410*4 = 1664) contains the seconds
 797    # we thus read four sectors continuously from the start of the file to
 798    # get maximum efficient I/O (remember we are trying to scan a lot of
 799    # files very quickly to build an index)
 800    fileobj.seek(0, os.SEEK_SET)
 801    block = fileobj.read(512 * 4)
 802
 803    # verify the signature with 4 bytes leading length
 804    signature = codecs.ascii_decode(block[4:12])[0]
 805    assert (signature == 'INTEHEAD')
 806
 807    # first record with the descriptor is 4 + 16 + 4 = 24 bytes long;
 808    # integer array is prefixed with another 4 bytes length, so the data
 809    # itself starts at 24 + 4 = 28 bytes (skipping 28 / 4 = 7 items --
 810    # they really didn't have aligned access in mind).
 811    # we don't care about the actual length of the block, as long as we
 812    # at least cover the number of items needed to extract the date
 813    intehead = numpy.frombuffer(block[28:],
 814                                dtype=numpy.dtype('>i4'),
 815                                count=512 - 7).copy()
 816
 817    # the data we read was big-endian; convert it into a native format
 818    # for this architecture
 819    if sys.byteorder != 'big':
 820        intehead.byteswap(True)
 821    intehead = intehead.view(numpy.int32)
 822
 823    # get the date for this header
 824    return _intehead_date(intehead)
 825
 826
 827class EclipseCase (object):
 828    """Read data for an Eclipse simulation case."""
 829
 830    def __init__(self, casename):
 831        """
 832        :param casename:  Path to the case, with or without extension.
 833        :type casename:   str
 834        """
 835        # keys of this dictionary will be dates, the values are the sequence
 836        # numbers that are associated with those dates
 837        self.by_date = {}
 838
 839        # get the directory of the file, using dot for current one if
 840        # no directory has been specified
 841        dir_name, file_name = os.path.split(casename)
 842        dir_name = '.' if dir_name == '' else path.expanduser(dir_name)
 843
 844        # further remove the extension from the filename
 845        root, _ = os.path.splitext(file_name)
 846
 847        # get a list of all restart files that are in this directory
 848        log.debug("Indexing restart files")
 849        for fname in os.listdir(dir_name):
 850            if fnmatch.fnmatch(fname,
 851                               '{0}.X[0-9][0-9][0-9][0-9]'.format(root)):
 852                # last four characters of the filename is the sequence number;
 853                # using int does not cause it to be interpreted octal even if
 854                # it starts with zero, fortunately
 855                seq = int(fname[-4:])
 856
 857                # quickly scan this file to determine the correlation between
 858                with open(os.path.join(dir_name, fname), 'rb') as fileobj:
 859                    this_date = _quick_date(fileobj)
 860                self.by_date[this_date] = seq
 861        log.debug("Found %d restart files", len(self.by_date))
 862
 863        # read the initial file to get the size of the grid and the
 864        # active flags, to load other properties
 865        self.root = os.path.join(dir_name, root)
 866        self.init = EclipseInit(self.root)
 867
 868        # we haven't loaded any grid or recurrent data yet
 869        self._grid = None
 870        self.recur = {}
 871        self.sum = {}
 872        self.comp = None
 873
 874    def shape(self):
 875        """\
 876        Get shape of returned field data.
 877
 878        :return:  (num_k, num_j, num_i)
 879        :rtype:   (int, int, int)
 880        """
 881        return self.init.shape
 882
 883    def start_date(self):
 884        """Starting date of the simulation"""
 885        return self.init.start_date
 886
 887    def components(self):
 888        """Components that exist in the restart file.
 889
 890        >>> case = EclipseCase (filename)
 891        >>> comps = case.components ()
 892        >>> print ("Number of components is %d" % (len (comps)))
 893        >>> for num, name in enumerate (comps):
 894        >>>     print ("%d : %s" % (num, name))
 895        """
 896        # get the last report date, and the index of that. the last
 897        # is selected because it is the most likely timestep to be
 898        # loaded by the user code (end results), and the number of
 899        # components shouldn't change during the run.
 900        lst_dat = max(self.by_date.keys())
 901        lst_seq = self.by_date[lst_dat]
 902        log.debug("Loading component list at date %s (step %04d)",
 903                  lst_dat.strftime('%Y-%m-%d'), lst_seq)
 904        restart = self._delay_load(lst_seq)
 905        return restart.components()
 906
 907    def phases(self):
 908        """Phases that exist in the restart file.
 909
 910        >>> case = EclipseCase (filename)
 911        >>> phs = case.phases ()
 912        >>> print ("Number of phases is %d" % (len (phs)))
 913        >>> print ("Keyword for first phase is %s" % ("S" + phs [0]))
 914        >>> if Phase.oil in phs:
 915        >>>     print ("Oil is present")
 916        >>> else:
 917        >>>     print ("Oil is not present")
 918        """
 919        return self.init.phases
 920
 921    def report_dates(self):
 922        """List of all the report dates that are available in this
 923        case. Items from this list can be used as a parameter to `at`
 924        to get the case for that particular report step.
 925
 926        :returns: List of available report dates in sequence
 927        :rtype:   list [datetime.datetime]
 928
 929        >>> for rstp in case.report_dates ():
 930                print (case.at (rstp).this_date)
 931
 932        .. seealso:: at
 933        """
 934        return sorted(self.by_date.keys())
 935
 936    def _delay_load(self, seq):
 937        """Make sure that recurrent data for a sequence step is loaded
 938        into memory.
 939
 940        :param seq:   Sequence number of the restart
 941        :type seq:    int
 942        """
 943        # check if we have loaded the data file yet. since we need
 944        # to index the file while reading from it, we maintain a
 945        # cache of all the files that has been loaded.
 946        if seq in self.recur:
 947            restart = self.recur[seq]
 948        else:
 949            restart = EclipseRestart(self.init, seq)
 950            self.recur[seq] = restart
 951
 952        return restart
 953
 954    def at(self, when):  # pylint: disable=invalid-name
 955        """Recurrent data for a certain timestep. The result of this
 956        method is usually passed to function that need to calculate
 957        something using several properties.
 958
 959        :param when:  Date of the property.
 960        :type when:   :class:`datetime.datetime` or int
 961        :returns:     Object containing properties for this timestep
 962        :rtype:       EclipseRestart
 963        """
 964        # get the sequence number of the report step; this is either
 965        # specified directly as an int, or given as a date
 966        if isinstance(when, datetime.datetime):
 967            assert (when in self.by_date)
 968            seq = self.by_date[when]
 969        else:
 970            seq = when
 971
 972        restart = self._delay_load(seq)
 973        return restart
 974
 975    def atsm(self, when):  # pylint: disable=invalid-name
 976        """Recurrent data from summary file for a certain timestep. The result
 977        of this method is usually passed to function that need to calculate
 978        something using several properties.
 979
 980        :param when:  Date of the property.
 981        :type when:   :class:`datetime.datetime` or int
 982        :returns:     Object containing summary properties for this timestep
 983        :rtype:       EclipseSummary
 984        """
 985        # get the sequence number of the report step; this is either
 986        # specified directly as an int, or given as a date
 987        if isinstance(when, datetime.datetime):
 988            assert (when in self.by_date)
 989            seq = self.by_date[when]
 990        else:
 991            seq = when
 992
 993        # check if we have loaded the data file yet. since we need
 994        # to index the file while reading from it, we maintain a
 995        # cache of all the files that has been loaded.
 996        if seq in self.sum:
 997            summary = self.sum[seq]
 998        else:
 999            summary = EclipseSummary(self.init, seq)
1000            self.sum[seq] = summary
1001
1002        return summary
1003
1004    def cell_data(self, prop, when=None):
1005        """Read cell-wise data from case. This can be either static information
1006        or recurrent information for a certain date.
1007
1008        :param prop:  Name of the property, e.g. 'SWAT'
1009        :type prop:   str
1010        :param when:  Date of the property, or None if static
1011        :type when:   :class:`datetime.datetime` or int
1012        :returns:     Loaded array for the property.
1013        :rtype:       :class:`numpy.ndarray`
1014
1015        :example:
1016        >>> case = EclipseCase (cmd_args.filename)
1017        >>> zmf2 = case.cell_data ('ZMF2', datetime.datetime (2054, 7, 1))
1018        """
1019        # if it is static property, it is found in the initial file
1020        if when is None:
1021            return self.init.cell_data(prop)
1022        else:
1023            # ask the restart file for the data
1024            return self.at(when).cell_data(prop)
1025
1026    def field_data(self, prop, when=None):
1027        """Read field-wise data from case. This can be either static information
1028        or recurrent information for a certain date.
1029
1030        :param prop:  Name of the property, e.g. 'ZPHASE'
1031        :type prop:   str
1032        :param when:  Date of the property, or None if static
1033        :type when:   :class:`datetime.datetime` or int
1034        :returns:     Loaded array for the property.
1035        :rtype:       :class:`numpy.ndarray`
1036
1037        :example:
1038        >>> case = EclipseCase (cmd_args.filename)
1039        >>> zphase = case.cell_data ('ZPHASE', datetime.datetime (2054, 7, 1))
1040        """
1041        # if it is static property, it is found in the initial file
1042        if when is None:
1043            return self.init.field_data(prop)
1044        else:
1045            # ask the restart file for the data
1046            return self.at(when).field_data(prop)
1047
1048    def summary_data(self, prop, when):
1049        """Read summary data from case. This is typically well data, but can
1050        also be e.g., newton iterations.
1051
1052        :param prop: Name of the property, e.g., 'WWPR PRO1'.
1053        :type prop:  str
1054        :param when: Date of the property
1055        :type when:  :class:`datetime.datetime` or int
1056        :returns:     Loaded array for the property.
1057        :rtype:       :class:`numpy.ndarray`
1058        :example:
1059        >>> case = EclipseCase (cmd_args.filename)
1060        >>> data = case.cell_data ('WWIR INJ-5',
1061                                   datetime.datetime (2054, 7, 1))
1062        """
1063        # Summary data is a bit more tricky than restart data. Mostly because
1064        # all specifications are given in the SMSPEC file. We must therefore
1065        # start by loading the mnemonic and the well associated with each data
1066        # vector
1067        return self.atsm(when).summary_data(prop)
1068
1069    def grid(self):
1070        """Grid structure for simulation case."""
1071        # on-demand load the grid file
1072        if self._grid is None:
1073            self._grid = EclipseGrid(self.init.root)
1074
1075        # offload this routine to the grid object
1076        return self._grid.grid()
1077
1078
1079class EclipseRFT (object):
1080    """Read data from an Eclipse RFT file."""
1081
1082    def __init__(self, casename):
1083        """
1084        :param casename:  Path to the case, with or without extension.
1085        :type casename:   str
1086        """
1087        # get the directory of the file, using dot for current one if
1088        # no directory has been specified
1089        dir_name, file_name = os.path.split(casename)
1090        dir_name = '.' if dir_name == '' else dir_name
1091
1092        # further remove the extension from the filename
1093        root, _ = os.path.splitext(file_name)
1094
1095        # read the initial file to get the size of the grid and the
1096        # active flags, to load other properties
1097        self.root = os.path.join(dir_name, root)
1098
1099        self.recur = {}
1100        self.sum = {}
1101
1102    def rft_data(self, well, prop):
1103        """
1104        Read the rft data for the requested well
1105        :param well: Name of well, e.g. 'PRO-1'
1106        :param prop: Type of property (depth, pressure, swat, or sgas)
1107        :type well: str
1108        :type prop: str
1109        :return: Loaded array for the property.
1110        :rtype: :class:`numpy.ndarray`
1111        >>> case = EclipseRFT (cmd_args.filename)
1112        >>> data = case.rft_data (well='INJ-5', prop='PRESSURE')
1113        """
1114        with EclipseFile(self.root, 'RFT') as eclf:
1115            vec_dat = None
1116            # Find the array from the current well
1117            for index in range(eclf.cnt['WELLETC ']):
1118                # well name for current index
1119                if eclf.get(kwd='WELLETC', seq=index)[1] == well.upper():
1120                    # check that this well has RFT data
1121                    assert 'R' in eclf.get(kwd='WELLETC', seq=index)[5]
1122                    vec_dat = eclf.get(kwd=prop, seq=index)
1123                    break
1124
1125        return vec_dat
1126
1127
1128def main(*args):
1129    """Read a data file to see if it parses OK."""
1130    # setup simple logging where we prefix each line with a letter code
1131    logging.basicConfig(level=logging.INFO,
1132                        format="%(levelname).1s: %(message).76s")
1133
1134    # parse command-line arguments
1135    parser = argparse.ArgumentParser()
1136    parser.add_argument("filename")
1137    cmd_args = parser.parse_args(*args)
1138
1139    # process file
1140    root, ext = os.path.splitext(cmd_args.filename)
1141    with EclipseFile(root, ext[1:]) as eclf:
1142        eclf.dump(True)
1143
1144
1145# executable library
1146if __name__ == "__main__":
1147    main(sys.argv[1:])