Enable subsets to be cut from a larger grdecl model

README.md

@@ -92,6 +92,8 @@ usage: em2ex.py [-h] [-o OUTPUT_FILE] [--filetype {eclipse,leapfrog}]
                 [--no-nodesets] [--no-sidesets] [-f] [-u] [--flip]
                 [--translate TRANSLATE TRANSLATE] [--mapaxes] [--pinch]
                 [--pinch-tol PINCH_TOL] [--refine-xy RX RY]
+                [--extract-i I_LO I_HI] [--extract-j J_LO J_HI]
+                [--extract-k K_LO K_HI]
                 filename
 
 Converts earth model to Exodus II format
@@ -121,6 +123,17 @@ options:
   --refine-xy RX RY     Refine the grid laterally by integer factors RX in x
                         and RY in y (vertical resolution unchanged). Each
                         child cell inherits its parent's element properties.
+  --extract-i I_LO I_HI
+                        Extract cells I_LO..I_HI along the x-axis (1-based
+                        inclusive, Eclipse-style). Cells are taken in file
+                        order, before any coordinate-system normalisation;
+                        runs before --refine-xy if both are given.
+  --extract-j J_LO J_HI
+                        Extract cells J_LO..J_HI along the y-axis (1-based
+                        inclusive).
+  --extract-k K_LO K_HI
+                        Extract cells K_LO..K_HI along the z-axis (1-based
+                        inclusive).
 ```
 
 ### Lateral refinement (Eclipse only)
@@ -133,6 +146,31 @@ The `--refine-xy RX RY` option refines an Eclipse grid in the (x, y) plane by in
 
 `RX` and `RY` must be strictly positive integers; anything else is rejected up front with an informative error.
 
+### Extracting a subset (Eclipse only)
+
+The `--extract-i`, `--extract-j` and `--extract-k` options pull a rectangular subset of cells out of a `grdecl` model along the x-, y- and z-axes respectively. Each takes two 1-based inclusive cell indices (Eclipse-style, matching the `BOX` keyword), and each is independently optional — any axis you don't restrict is kept in full. For example, to keep only cells `i=10..30, j=5..40` across every layer:
+
+```bash
+./em2ex.py --extract-i 10 30 --extract-j 5 40 large.grdecl
+```
+
+Or to take only the top 20 layers:
+
+```bash
+./em2ex.py --extract-k 1 20 large.grdecl
+```
+
+Indices refer to cell positions **as they appear in the file** — same numbering you see in the `SPECGRID` keyword and in the order properties like `PORO` are listed. The subset is taken before any further coordinate processing.
+
+**Composition with `--refine-xy`.** If both are given, extract runs first and refinement applies to the subset (so `--extract-i 10 30 --refine-xy 2 2` extracts 21 cells along the x-axis and then refines to 42; it does *not* refine the full grid and then extract from it). This is almost always what you want — refining the whole grid just to throw most of it away would be wasteful.
+
+**Composition with `--flip` and with left-hand coordinate files.** Extract operates on the file's i/j/k indexing *before* `em2ex` flips the z values (`--flip`) or normalises a left-handed coordinate system to right-handed (which happens automatically when the file's x or y coordinates decrease). The practical consequence:
+
+- `--flip` doesn't affect which cells are extracted — only their z sign in the output. The k range you give is the same range you'd give without `--flip`.
+- For a left-handed coordinate file, the extracted region is still the cells at file indices `i=I_LO..I_HI` (etc.), exactly as they're listed in `PORO`. In the output mesh, those cells then get re-numbered through the auto-flip to the canonical right-handed system, so their i (and/or j) indices in the produced Exodus mesh may run in the opposite direction from the file's. The geometry and properties are preserved; only the index sense changes.
+
+If any range is out of bounds for the file's `SPECGRID` size, or if `LO > HI`, the conversion is rejected up front with the actual dimensions cited.
+
 `em2ex` attempts to guess the reservoir model format from the file extension (see supported formats below). If the reservoir model has a non-standard file extension, the user can force
 `em2ex` to read the correct format using the `--filetype` commandline option.
 

em2ex.py

@@ -9,15 +9,16 @@
 import os
 
 def _positive_int(s):
-    ''' argparse type for strictly positive integer refinement factors. '''
+    ''' argparse type for strictly positive integers (refinement factors and
+    extract indices). '''
     try:
         v = int(s)
     except (TypeError, ValueError):
         raise argparse.ArgumentTypeError(
-            "refinement factor must be a positive integer, got {!r}".format(s))
+            "expected a positive integer, got {!r}".format(s))
     if v < 1:
         raise argparse.ArgumentTypeError(
-            "refinement factor must be >= 1, got {}".format(v))
+            "expected a positive integer (>= 1), got {}".format(v))
     return v
 
 def get_parser():
@@ -40,6 +41,15 @@ def get_parser():
     parser.add_argument('--refine-xy', nargs = 2, dest = 'refine_xy', type = _positive_int,
         metavar = ('RX', 'RY'),
         help = 'Refine the grid laterally by integer factors RX in x and RY in y (vertical resolution unchanged). Each child cell inherits its parent\'s element properties.')
+    parser.add_argument('--extract-i', nargs = 2, dest = 'extract_i', type = _positive_int,
+        metavar = ('I_LO', 'I_HI'),
+        help = 'Extract cells I_LO..I_HI along the x-axis (1-based inclusive, Eclipse-style). Cells are taken in file order, before any coordinate-system normalisation; runs before --refine-xy if both are given.')
+    parser.add_argument('--extract-j', nargs = 2, dest = 'extract_j', type = _positive_int,
+        metavar = ('J_LO', 'J_HI'),
+        help = 'Extract cells J_LO..J_HI along the y-axis (1-based inclusive).')
+    parser.add_argument('--extract-k', nargs = 2, dest = 'extract_k', type = _positive_int,
+        metavar = ('K_LO', 'K_HI'),
+        help = 'Extract cells K_LO..K_HI along the z-axis (1-based inclusive).')
     return parser
 
 def main():

readers/eclipse.py

@@ -240,6 +240,27 @@ def parseEclipse(f, args):
     # The COORD data has six entries for each of the (nx+1)*(ny+1) nodes
     coord = np.asarray(eclipse.coord).reshape(ny+1, nx+1, 6)
 
+    # Reshape ZCORN early so subsetting (--extract-*) can slice it. The reshape
+    # is independent of any later coord transforms (flip / translate / mapaxes /
+    # flip_z), all of which leave the (k, j, i) cell indexing unchanged.
+    zcorn = np.asarray(eclipse.zcorn).reshape(2*nz, 2*ny, 2*nx)
+
+    # Apply --extract-i/-j/-k subsetting if requested. Indices are 1-based
+    # inclusive in file order (matching the cells as they appear in the
+    # grdecl SPECGRID / properties sections), and the slice happens before
+    # any flip / translate / mapaxes / refine so the user never has to think
+    # about coordinate-system normalisation. See README for the flip caveat.
+    extract_i = getattr(args, 'extract_i', None)
+    extract_j = getattr(args, 'extract_j', None)
+    extract_k = getattr(args, 'extract_k', None)
+    if extract_i or extract_j or extract_k:
+        i_lo, i_hi = _resolve_extract_range(extract_i, nx, 'i')
+        j_lo, j_hi = _resolve_extract_range(extract_j, ny, 'j')
+        k_lo, k_hi = _resolve_extract_range(extract_k, nz, 'k')
+        coord, zcorn, eclipse.elemProps, nx, ny, nz = extractSubgrid(
+            coord, zcorn, eclipse.elemProps, nx, ny, nz,
+            i_lo, i_hi, j_lo, j_hi, k_lo, k_hi)
+
     # The exodus node numbering relies on a right-hand coordinate system, with
     # x and y increasing. However, eclipse can output a grid with a left-hand
     # coordinate system, with either (or both) x and y decreasing (ie, pointing in
@@ -293,11 +314,6 @@ def parseEclipse(f, args):
                 coord[:,:,xi] = xorigin + xdata * xvec[0] + ydata * yvec[0]
                 coord[:,:,yi] = yorigin + xdata * xvec[1] + ydata * yvec[1]
 
-    # The ZCORN data varies by x, y and then z. Reshape it into an array where
-    # the first index gives the layer, the second the y coordinates and the third
-    # the x coordinates
-    zcorn = np.asarray(eclipse.zcorn).reshape(2*nz, 2*ny, 2*nx)
-
     # Flip the Z coordinates if specified
     if args.flip_z:
         zcorn = - zcorn
@@ -595,6 +611,54 @@ def numberNodesInElems(elemcornz, active_elements):
     return elemNodes_flat.reshape(nz, ny, nx, 8).astype(int)
 
 
+def _resolve_extract_range(rng, n, range_letter):
+    ''' Validate a 1-based inclusive extract range against the file's axis
+    size and return the 0-based half-open [lo, hi) pair to slice with. None
+    means "keep the full axis". `range_letter` is the lowercase i/j/k that
+    appears in the CLI flag and the SPECGRID NX/NY/NZ symbol. '''
+    if rng is None:
+        return 0, n
+    axis = {'i': 'x', 'j': 'y', 'k': 'z'}[range_letter]
+    specgrid_dim = {'i': 'NX', 'j': 'NY', 'k': 'NZ'}[range_letter]
+    lo_1, hi_1 = rng
+    if lo_1 > hi_1:
+        print("--extract-{} requires LO <= HI, got {} > {}".format(
+            range_letter, lo_1, hi_1))
+        exit()
+    if lo_1 < 1 or hi_1 > n:
+        print("--extract-{} range {}..{} is out of bounds for the {}-axis ({}={} in SPECGRID)".format(
+            range_letter, lo_1, hi_1, axis, specgrid_dim, n))
+        exit()
+    return lo_1 - 1, hi_1
+
+
+def extractSubgrid(coord, zcorn, elemProps, nx, ny, nz,
+                   i_lo, i_hi, j_lo, j_hi, k_lo, k_hi):
+    ''' Slice the grid to the given (i, j, k) ranges. Inputs are in file order
+    (the slice runs before any flip / translate / mapaxes); ranges are 0-based
+    half-open [lo, hi).
+
+    Inputs:
+        coord     : (ny+1, nx+1, 6)        pillar top/bottom (x, y, z)
+        zcorn     : (2*nz, 2*ny, 2*nx)     per-cell corner z values
+        elemProps : dict[name -> flat ndarray of length nx*ny*nz]
+
+    Returns (coord, zcorn, elemProps, new_nx, new_ny, new_nz) for the subset.
+    '''
+    # Pillars in (j, i): keep one extra pillar past the high end to bound the
+    # last cell. .copy() decouples from the original so downstream writes (e.g.
+    # the translate block) don't accidentally mutate it.
+    coord = coord[j_lo:j_hi+1, i_lo:i_hi+1, :].copy()
+    zcorn = zcorn[2*k_lo:2*k_hi, 2*j_lo:2*j_hi, 2*i_lo:2*i_hi].copy()
+
+    new_elemProps = {}
+    for name, vals in elemProps.items():
+        p = np.asarray(vals).reshape(nz, ny, nx)[k_lo:k_hi, j_lo:j_hi, i_lo:i_hi]
+        new_elemProps[name] = p.flatten()
+
+    return coord, zcorn, new_elemProps, i_hi - i_lo, j_hi - j_lo, k_hi - k_lo
+
+
 def refineLaterally(coord, zcorn, elemProps, nx, ny, nz, rx, ry):
     ''' Refine the grid laterally by integer factors (rx, ry) without refining
     vertically. Pillars (and their endpoint x, y, z stored in coord) are linearly

test/eclipse/simple_cube_extract.grdecl

@@ -0,0 +1,80 @@
+-- 4x3x2 grid with uniquely-valued PORO per cell.
+-- Combined with --extract-i 2 3 --extract-j 1 2 --extract-k 1 1 this should
+-- produce a 2x2x1 grid geometrically and per-property equivalent to
+-- simple_cube_extract_ref.grdecl. PORO values are 0.<k><j><i> (1-based) so
+-- you can read which cell survived the extract directly off the values.
+
+SPECGRID
+4 3 2 1 F /
+
+GRIDUNIT
+  METRES /
+
+COORD
+ 0.000 0.000 0.000 0.000 0.000 2.000
+ 1.000 0.000 0.000 1.000 0.000 2.000
+ 2.000 0.000 0.000 2.000 0.000 2.000
+ 3.000 0.000 0.000 3.000 0.000 2.000
+ 4.000 0.000 0.000 4.000 0.000 2.000
+ 0.000 1.000 0.000 0.000 1.000 2.000
+ 1.000 1.000 0.000 1.000 1.000 2.000
+ 2.000 1.000 0.000 2.000 1.000 2.000
+ 3.000 1.000 0.000 3.000 1.000 2.000
+ 4.000 1.000 0.000 4.000 1.000 2.000
+ 0.000 2.000 0.000 0.000 2.000 2.000
+ 1.000 2.000 0.000 1.000 2.000 2.000
+ 2.000 2.000 0.000 2.000 2.000 2.000
+ 3.000 2.000 0.000 3.000 2.000 2.000
+ 4.000 2.000 0.000 4.000 2.000 2.000
+ 0.000 3.000 0.000 0.000 3.000 2.000
+ 1.000 3.000 0.000 1.000 3.000 2.000
+ 2.000 3.000 0.000 2.000 3.000 2.000
+ 3.000 3.000 0.000 3.000 3.000 2.000
+ 4.000 3.000 0.000 4.000 3.000 2.000
+/
+
+ZCORN
+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 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 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 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
+2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
+2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
+2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
+2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
+2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
+/
+
+ACTNUM
+1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1
+/
+
+PORO
+0.111 0.112 0.113 0.114
+0.121 0.122 0.123 0.124
+0.131 0.132 0.133 0.134
+0.211 0.212 0.213 0.214
+0.221 0.222 0.223 0.224
+0.231 0.232 0.233 0.234
+/
+
+SATNUM
+1 1 1 1 1 1 1 1 1 1 1 1
+1 1 1 1 1 1 1 1 1 1 1 1
+/

test/eclipse/tests

@@ -68,6 +68,12 @@ simple_cube_refine:
   cli_args: --refine-xy 2 3
   gold: simple_cube_refine.e
 
+simple_cube_extract:
+  filename: simple_cube_extract.grdecl
+  type: exodiff
+  cli_args: --extract-i 2 3 --extract-j 1 2 --extract-k 1 1
+  gold: simple_cube_extract.e
+
 missing_specgrid:
   filename: missing_specgrid.grdecl
   type: exception
@@ -102,22 +108,46 @@ refine_xy_zero:
   filename: simple_cube.grdecl
   type: exception
   cli_args: --refine-xy 0 2
-  expected_error: refinement factor must be >= 1
+  expected_error: expected a positive integer (>= 1)
 
 refine_xy_negative:
   filename: simple_cube.grdecl
   type: exception
   cli_args: --refine-xy 2 -3
-  expected_error: refinement factor must be >= 1
+  expected_error: expected a positive integer (>= 1)
 
 refine_xy_non_integer:
   filename: simple_cube.grdecl
   type: exception
   cli_args: --refine-xy 1.5 2
-  expected_error: refinement factor must be a positive integer
+  expected_error: expected a positive integer
 
 refine_xy_garbage:
   filename: simple_cube.grdecl
   type: exception
   cli_args: --refine-xy foo 2
-  expected_error: refinement factor must be a positive integer
+  expected_error: expected a positive integer
+
+extract_i_out_of_bounds:
+  filename: simple_cube.grdecl
+  type: exception
+  cli_args: --extract-i 1 99
+  expected_error: --extract-i range 1..99 is out of bounds for the x-axis (NX=3 in SPECGRID)
+
+extract_j_inverted:
+  filename: simple_cube.grdecl
+  type: exception
+  cli_args: --extract-j 3 1
+  expected_error: --extract-j requires LO <= HI, got 3 > 1
+
+extract_k_zero:
+  filename: simple_cube.grdecl
+  type: exception
+  cli_args: --extract-k 0 2
+  expected_error: expected a positive integer (>= 1)
+
+extract_i_non_integer:
+  filename: simple_cube.grdecl
+  type: exception
+  cli_args: --extract-i 1.5 2
+  expected_error: expected a positive integer