misc.grid.sector

  1"""\
  2Extract a sector from an existing cornerpoint grid.
  3"""
  4import argparse
  5import collections
  6import logging
  7import numpy
  8import re
  9import sys
 10
 11import misc.grid as pyr
 12import misc.grdecl as grdecl
 13
 14
 15# add a valid log in case we are not run through the main program which
 16# sets up one for us
 17log = logging.getLogger(__name__)  # pylint: disable=invalid-name
 18log.addHandler(logging.NullHandler())
 19
 20# three consecutive integers, separated by comma, and perhaps some spaces
 21# thrown in for readability, optionally enclosed in parenthesis
 22tuple_format = re.compile(r'\(?([0-9]+)\ *\,\ *([0-9]+)\ *\,\ *([0-9]+)\)?')
 23
 24
 25def parse_tuple(corner):
 26    """\
 27    Parse a parenthesized tuple into three constituent coordinates.
 28
 29    :param corner:  Coordinate specification on the format "(i1,j1,k1)"
 30    :type  corner:  str
 31    :returns:       The parsed tuple, converted into zero-based coordinates
 32                    and in Python-matrix order: (k, j, i)
 33    :rtype  :       (int, int, int)
 34    """
 35    # let the regular expression engine parse the string
 36    match = re.match(tuple_format, corner.strip())
 37
 38    # if the string matched, then we know that each of the sub-groups can be
 39    # parsed into strings successfully. group 0 is the entire string, so we
 40    # get natural numbering into the parenthesized expressions. subtract one
 41    # to get zero-based coordinates
 42    if match:
 43        i = int(match.group(1)) - 1
 44        j = int(match.group(2)) - 1
 45        k = int(match.group(3)) - 1
 46        return (k, j, i)
 47    # if we didn't got any valid string, then return a bottom value
 48    else:
 49        return None
 50
 51
 52def sort_tuples(corner, opposite):
 53    """\
 54    :param corner:    Coordinates of one corner
 55    :type  corner:    (int, int, int)
 56    :param opposite:  Coordinates of the opposite corner
 57    :type  opposite:  (int, int, int)
 58    :returns:         The two tuples, but with coordinates interchanged so that
 59                      one corner is always in the lower, left, back and the
 60                      other is in the upper, right, front
 61    :rtype:           ((int, int, int), (int, int, int))
 62    """
 63    # pick out the most extreme variant in either direction, into each its own
 64    # variable; this may be the same as the input or not, but at least we know
 65    # for sure when we return from this method
 66    least = (min(corner[0], opposite[0]),
 67             min(corner[1], opposite[1]),
 68             min(corner[2], opposite[2]))
 69    most = (max(corner[0], opposite[0]),
 70            max(corner[1], opposite[1]),
 71            max(corner[2], opposite[2]))
 72    return (least, most)
 73
 74
 75def extract_dimens(least, most):
 76    """\
 77    Build a new dimension tuple for a submodel
 78
 79    :param least:  Lower, left-most, back corner of submodel, (k1, j1, i1)
 80    :type  least:  (int, int, int)
 81    :param most:   Upper, right-most, front corner of submodel, (k2, j2, i2)
 82    :type  most:   (int, int, int)
 83    :returns:      Dimensions of the submodel
 84    :rtype:        numpy.ndarray((3,))
 85    """
 86    # split the corners into constituents
 87    k1, j1, i1 = least
 88    k2, j2, i2 = most
 89
 90    # make an array out of the cartesian distance of the two corners
 91    sector_dimens = numpy.array([i2-i1+1, j2-j1+1, k2-k1+1], dtype=numpy.int32)
 92    return sector_dimens
 93
 94
 95def extract_coord(coord, least, most):
 96    """\
 97    Extract the coordinate pillars for a submodel
 98
 99    :param coord:  Coordinate pillars for the entire grid
100    :type  coord:  numpy.ndarray((nj+1, ni+1, 2, 3))
101    :param least:  Lower, left-most, back corner of submodel, (k1, j1, i1)
102    :type  least:  (int, int, int)
103    :param most:   Upper, right-most, front corner of submodel, (k2, j2, i2)
104    :type  most:   (int, int, int)
105    :returns:      Coordinate pilars for the submodel
106    :rtype:        numpy.ndarray((j2-j1+2, i2-i1+2, 2, 3))
107    """
108    # split the corners into constituents
109    k1, j1, i1 = least
110    k2, j2, i2 = most
111
112    # add one to get the pillar on the other side of the element, so that
113    # we include the last element, add one more since Python indexing is
114    # end-exclusive
115    sector_coord = coord[j1:(j2+2), i1:(i2+2), :, :]
116    return sector_coord
117
118
119def extract_zcorn(zcorn, least, most):
120    """\
121    Extract the hinge depth values for a submodel
122
123    :param zcorn:  Hinge depth values for the entire grid
124    :type  zcorn:  numpy.ndarray((nk, 2, nj, 2, ni, 2))
125    :param least:  Lower, left-most, back corner of submodel, (k1, j1, i1)
126    :type  least:  (int, int, int)
127    :param most:   Upper, right-most, front corner of submodel, (k2, j2, i2)
128    :type  most:   (int, int, int)
129    :returns:      Hinge depth values for the submodel
130    :rtype:        numpy.ndarray((k2-k1+1, 2, j2-j1+1, 2, i2-i1+1))
131    """
132    # split the corners into constituents
133    k1, j1, i1 = least
134    k2, j2, i2 = most
135
136    # add one since Python-indexing is end-exclusive, and we want to include
137    # the element in the opposite corner. all eight hinges are returned for
138    # each element (we only select in every other dimension)
139    sector_zcorn = zcorn[k1:(k2+1), :, j1:(j2+1), :, i1:(i2+1), :]
140    return sector_zcorn
141
142
143def extract_cell_prop(prop, least, most):
144    """\
145    Extract the property values for a submodel
146
147    :param prop:   Property values for each cell in the the entire grid
148    :type  prop:   numpy.ndarray((nk, nj, ni))
149    :param least:  Lower, left-most, back corner of submodel, (k1, j1, i1)
150    :type  least:  (int, int, int)
151    :param most:   Upper, right-most, front corner of submodel, (k2, j2, i2)
152    :type  most:   (int, int, int)
153    :returns:      Property values for each cell in the submodel
154    :rtype:        numpy.ndarray((k2-k1+1, j2-j1+1, i2-i1+1))
155    """
156    # split the corners into constituents
157    k1, j1, i1 = least
158    k2, j2, i2 = most
159
160    # add one since Python-indexing is end-exclusive, and we want to include
161    # the element in the opposite corner.
162    sector_prop = prop[k1:(k2+1), j1:(j2+1), i1:(i2+1)]
163    return sector_prop
164
165
166def extract_grid(grid, least, most):
167    """\
168    Extract a submodel from a full grid
169
170    :param grid:   Attributes of the full grid, like COORD, ZCORN, ACTNUM
171    :type  grid:   dict
172    :param least:  Lower, left-most, back corner of submodel, (k1, j1, i1)
173    :type  least:  (int, int, int)
174    :param most:   Upper, right-most, front corner of submodel, (k2, j2, i2)
175    :type  most:   (int, int, int)
176    :returns:      Attributes of the sector model
177    :rtype:        dict
178    """
179    # create a new grid and fill extract standard properties
180    sector = collections.OrderedDict()
181    sector['DIMENS'] = extract_dimens(least, most)
182    sector['COORD'] = extract_coord(grid['COORD'], least, most)
183    sector['ZCORN'] = extract_zcorn(grid['ZCORN'], least, most)
184    sector['ACTNUM'] = extract_cell_prop(grid['ACTNUM'], least, most)
185
186    # then extract all extra cell properties, such as PORO, PERMX that
187    # may have been added
188    for prop_name in grid:
189        # ignore the standard properties; they have already been converted
190        # by specialized methods above
191        prop_upper = prop_name.upper()
192        std_prop = ((prop_upper == 'DIMENS') or (prop_upper == 'COORD') or
193                    (prop_upper == 'ZCORN') or (prop_upper == 'ACTNUM'))
194        # use the generic method to convert this property
195        if not std_prop:
196            sector[prop_name] = extract_cell_prop(grid[prop_name], least, most)
197
198    return sector
199
200
201def main(*args):
202    """Read a data file to see if it parses OK."""
203    # setup simple logging where we prefix each line with a letter code
204    logging.basicConfig(level=logging.INFO,
205                        format="%(levelname).1s: %(message).76s")
206
207    # parse command-line arguments
208    parser = argparse.ArgumentParser()
209    parser.add_argument("infile", metavar="infile.grdecl", help="Input model")
210    parser.add_argument("outfile", metavar="outfile",
211                        help="Output sector model")
212    parser.add_argument("first", metavar="i1,j1,k1", type=parse_tuple,
213                        help="A corner of the sector model (one-based)")
214    parser.add_argument("last", metavar="i2,j2,k2", type=parse_tuple,
215                        help="The opposite corner of the sector (one-based)")
216    parser.add_argument("--dialect", choices=['ecl', 'cmg'], default='ecl')
217    parser.add_argument("--verbose", action='store_true')
218    parser.add_argument("--quiet", action='store_true')
219    cmd_args = parser.parse_args(*args)
220
221    # adjust the verbosity of the program
222    if cmd_args.verbose:
223        logging.getLogger(__name__).setLevel(logging.DEBUG)
224    if cmd_args.quiet:
225        logging.getLogger(__name__).setLevel(logging.NOTSET)
226
227    # read the input file
228    log.info('Reading full grid')
229    grid = pyr.read_grid(cmd_args.infile)
230
231    # get the two opposite corners that defines the submodel
232    log.info('Determining scope of sector model')
233    least, most = sort_tuples(cmd_args.first, cmd_args.last)
234    log.info('Sector model is (%d, %d, %d)-(%d, %d, %d)',
235             least[2]+1, least[1]+1, least[0]+1,
236             most[2]+1, most[1]+1, most[0]+1)
237
238    # extract the data for the submodel into a new grid
239    log.info('Extracting sector model from full grid')
240    submodel = extract_grid(grid, least, most)
241
242    # write this grid to output
243    log.info('Writing sector model to disk')
244    grdecl.write(cmd_args.outfile, submodel, cmd_args.dialect,
245                 multi_file=True)
246
247
248# if this module is called as a standalone program, then pass all the
249# arguments, except the program name
250if __name__ == "__main__":
251    main(sys.argv[1:])