misc.grid.cornerpoint

  1"""Common functionality for visualization of corner-point grids.
  2
  3import pyresito.grid.cornerpoint as cp
  4"""
  5import itertools as it
  6import logging
  7import numpy as np
  8import numpy.ma as ma
  9
 10
 11# add a valid log in case we are not run through the main program which
 12# sets up one for us
 13log = logging.getLogger(__name__)  # pylint: disable=invalid-name
 14log.addHandler(logging.NullHandler())
 15
 16
 17def scatter(cell_field):
 18    """\
 19    Duplicate all items in every dimension.
 20
 21    Use this method to duplicate values associated with each cell to
 22    each of the corners in the cell.
 23
 24    :param cell_field:   Property per cell in the grid
 25    :type  cell_field:   :class:`numpy.ndarray`, shape = (nk, nj, ni)
 26    :return:             Property per corner in the cube
 27    :rtype:              :class:`numpy.ndarray`, shape = (nk, 2, nj, 2, ni, 2)
 28
 29    >>> scatter(np.array([[[1, 2],
 30                           [3, 4]],
 31                          [[5, 6],
 32                           [7, 8]]]))
 33
 34    array([[[[[[1, 1], [2, 2]],
 35              [[1, 1], [2, 2]]],
 36             [[[3, 3], [4, 4]],
 37              [[3, 3], [4, 4]]]],
 38
 39            [[[[1, 1], [2, 2]],
 40              [[1, 1], [2, 2]]],
 41             [[[3, 3], [4, 4]],
 42              [[3, 3], [4, 4]]]]],
 43
 44           [[[[[5, 5], [6, 6]],
 45              [[5, 5], [6, 6]]],
 46             [[[7, 7], [8, 8]],
 47              [[7, 7], [8, 8]]]],
 48
 49            [[[[5, 5], [6, 6]],
 50              [[5, 5], [6, 6]]],
 51             [[[7, 7], [8, 8]],
 52              [[7, 7], [8, 8]]]]]])
 53    """
 54    # get the dimensions of the cube in easily understood aliases
 55    nk, nj, ni = cell_field.shape  # pylint: disable=invalid-name
 56
 57    # create an extra dimension so that we can duplicate individual values
 58    flat = np.reshape(cell_field, (nk, nj, ni, 1))
 59
 60    # duplicate along each of the axis
 61    dup_i = np.tile(flat, (1, 1, 1, 2))
 62    dup_j = np.tile(dup_i, (1, 1, 2, 1))
 63    dup_k = np.tile(dup_j, (1, 2, 1, 1))
 64
 65    # reformat to a cube with the appropriate dimensions
 66    corn_field = np.reshape(dup_k, (nk, 2, nj, 2, ni, 2))
 67    return corn_field
 68
 69
 70def inner_dup(pillar_field):
 71    """Duplicate all inner items in both dimensions.
 72
 73    Use this method to duplicate values associated with each pillar to
 74    the corners on each side(s) of the pillar; four corners for all
 75    interior pillars, two corners for all pillars on the rim and only
 76    one (element) corner for those pillars that are on the grid corners.
 77
 78    :param pillar_field: Property per pillar in the grid
 79    :type  pillar_field: :class:`numpy.ndarray`, shape=(m+1, n+1)
 80    :returns:            Property per corner in a grid plane
 81    :rtype:              :class:`numpy.ndarray`, shape=(2*m, 2*n))
 82
 83    >>> inner_dup(np.array ([[1, 2, 3],
 84                             [4, 5, 6],
 85                             [7, 8, 9]]))
 86    array([[1, 2, 2, 3],
 87           [4, 5, 5, 6],
 88           [4, 5, 5, 6],
 89           [7, 8, 8, 9]])
 90    """
 91    # extract the inner horizontal part of the plane
 92    horz_part = pillar_field[:, 1:-1]
 93
 94    # tile it to a separate plane
 95    horz_dupd = np.tile(horz_part, (2, 1, 1))
 96
 97    # shift the tile to the end so it can be rolled
 98    # into the second dimension
 99    horz_shft = np.transpose(horz_dupd, (1, 2, 0))
100
101    # roll this into the second dimension
102    horz_roll = np.reshape(horz_shft,
103                           (horz_shft.shape[0],
104                            horz_shft.shape[1] * 2))
105
106    # add back the first and last column
107    horz = np.column_stack([pillar_field[:, 0],
108                            horz_roll,
109                            pillar_field[:, -1]])
110
111    # extract the inner vertical part of the plane
112    vert_part = horz[1:-1, :]
113
114    # tile it to a separate plane
115    vert_dupd = np.tile(vert_part, (1, 1, 2))
116
117    # roll this into the first dimension
118    vert_roll = np.reshape(vert_dupd,
119                           (vert_dupd.shape[1] * 2,
120                            vert_dupd.shape[2] // 2))
121
122    # add back the first and last row
123    result = np.vstack([horz[0, :],
124                        vert_roll,
125                        horz[-1, :]])
126    return result
127
128
129def elem_vtcs_ndcs(nk, nj, ni):  # pylint: disable=invalid-name
130    """List zcorn indices used by every element in an nk*nj*ni grid
131
132    :param nk: Number of layers in Z direction
133    :type nk:  int
134    :param nj: Number of elements in the Y direction
135    :type nj:  int
136    :param ni: Number of elements in the X direction
137    :type ni:  int
138    :returns:  Zero-based indices for the hexahedral element corners
139    :rtype:    :class:`numpy.ndarray`, shape=(nk*nj*ni, 8), dtype=int
140    """
141    # hex_perm is the order a hexahedron should be specified to the
142    # gridding engine (it is the same for VTK and Ensight Gold formats)
143    # relative to the order they are in the Eclipse format. the latter
144    # four is listed first because we turn the z order around with a
145    # transform.
146    # local indices:      0, 1, 2, 3, 4, 5, 6, 7
147    hex_perm = np.array([4, 5, 7, 6, 0, 1, 3, 2], dtype=int)
148
149    # kji is the global (k, j, i) index for every element; we tile
150    # this into 8 corners per element, with 3 indices per corner;
151    # this this matrix contains repeats of 8 identical permutations
152    # of the values 0..(ni-1), 0..(nj-1) and 0..(nk-1)
153    kji = np.array(list(it.product(range(nk),
154                                   range(nj),
155                                   range(ni))))
156    kji = np.reshape(np.tile(kji, [1, hex_perm.shape[0]]), [-1, 3])
157
158    # bfr is the local (bottom, front, right) index for every corner;
159    # this is a block of 8 local corners which are picked out of the
160    # zcorn matrix to then correspond with the gridding engines order.
161    # we tile this block for every element
162    bfr = np.array(list(it.product([0, 1], repeat=3)))[hex_perm, :]
163    bfr = np.tile(bfr, [nk*nj*ni, 1])
164
165    # calculate running number for every corner, combining the element
166    # indices kji with the local corner indices bfr. the zcorn hypercube
167    # has dimensions nk * 2 * nj * 2 * ni * 2; the indexing is simply a
168    # flattened view of this hypercube.
169    ndx = (((2*kji[:, 0] + bfr[:, 0]) * 2*nj +
170            (2*kji[:, 1] + bfr[:, 1])) * 2*ni +
171           2*kji[:, 2] + bfr[:, 2])
172
173    # we then collect the indices so that we get 8 of them (all those that
174    # belongs to the same element) on a row
175    ndx = np.reshape(ndx, [-1, hex_perm.shape[0]])
176    return ndx
177
178
179# pylint: disable=invalid-name, multiple-statements, too-many-locals
180def corner_coordinates(coord, zcorn):
181    """Generate (x, y, z) coordinates for each corner-point.
182
183    :param coord: Pillar geometrical information
184    :type coord:  :class:`numpy.ndarray`, shape=(nj+1, ni+1, 2, 3)
185    :param zcorn: Depth values along each pillar
186    :type zcorn:  :class:`numpy.ndarray`, shape=(nk, 2, nj, 2, ni, 2)
187    :returns:     Coordinate values for each corner
188    :rtype:       :class:`numpy.ndarray`, shape=(3, nk*2*nj*2*ni*2)
189    """
190    # indices to treat the pillar matrix as a record
191    X = 0
192    Y = 1
193    Z = 2
194    START = 0
195    END = 1
196
197    # calculate depth derivatives for the pillar slope
198    dx = coord[:, :, END, X] - coord[:, :, START, X]
199    dy = coord[:, :, END, Y] - coord[:, :, START, Y]
200    dz = coord[:, :, END, Z] - coord[:, :, START, Z]
201
202    # make the pillar information compatible with a plane of corners,
203    # having a singleton dimension at the front (instead tiling 2*nk)
204    x0 = inner_dup(coord[:, :, START, X])[None, :, :]
205    y0 = inner_dup(coord[:, :, START, Y])[None, :, :]
206    z0 = inner_dup(coord[:, :, START, Z])[None, :, :]
207    dxdz = inner_dup(dx / dz)[None, :, :]
208    dydz = inner_dup(dy / dz)[None, :, :]
209
210    # infer the grid size from the dimensions of coord
211    nj = coord.shape[0] - 1
212    ni = coord.shape[1] - 1
213
214    # make the formal dimensions of the depths compatible with the plane
215    # of vectors created from the pillars; notice that the first dimension
216    # is not singular
217    z = np.reshape(zcorn, (-1, nj*2, ni*2))
218
219    # calculate the point in the plane where the horizon intersect with the
220    # pillar; we get the (x, y)-coordinates for each corner. notice that the
221    # difference z - z0 is for each point, whereas dz is for the pillar!
222    x = (z - z0) * dxdz + x0
223    y = (z - z0) * dydz + y0
224
225    # reformat the coordinates into the final coordinates
226    xyz = np.vstack([np.reshape(x, (1, -1)),
227                     np.reshape(y, (1, -1)),
228                     np.reshape(z, (1, -1))])
229    return xyz
230
231
232# Enumeration of the dimensions in a point tuple
233Dim = type('Dim', (), {
234    'X': 0,
235    'Y': 1,
236    'Z': 2,
237})
238
239
240# Enumeration of the faces of a hexahedron. The contents are the local
241# indices in each row of 8 corners that are set up for the elements by
242# the elem_vtcs_ndcs function. (The indices here are the numbers that
243# are in the comment "local indices" above the hex_perm selector)
244Face = type('Face', (), {
245    'DOWN':  np.array([0, 1, 2, 3]),
246    'UP':    np.array([4, 5, 6, 7]),
247    'FRONT': np.array([4, 7, 3, 0]),
248    'BACK':  np.array([5, 6, 2, 1]),
249    'LEFT':  np.array([7, 6, 2, 3]),
250    'RIGHT': np.array([4, 5, 1, 0]),
251    'ALL':   np.array([0, 1, 2, 3, 4, 5, 6, 7]),
252})
253
254
255def cp_cells(grid, face):
256    """
257    Make a cell array from a cornerpoint grid. The cells will be in the
258    same order as the Cartesian enumeration of the grid.
259
260    :param grid: Pillar coordinates and corner depths
261    :type grid:  dict, with 'coord' and 'zcorn' properties
262    :param face: Which face that should be extracted, e.g. Face.UP
263    :type face:  Face enum
264    :returns:    Set of geometrical objects that can be sent to rendering
265    :rtype:      dict with 'points', shape (nverts, 3), and 'cells', shape
266                 (nelems, ncorns), where ncorns is either 8 (hexahedron
267                 volume or 4 (quadrilateral face), depending on the face
268                 parameter.
269    """
270    src = {}
271
272    log.debug("Generating points for corner depths")
273    # notice that VTK wants shape (nverts, 3), where as the code sets up
274    # for Ensight format which is (3, nverts).
275    src['points'] = corner_coordinates(grid['COORD'],
276                                       grid['ZCORN']).T.copy()
277
278    # order of the dimensions are different from the structure to the call
279    ni, nj, nk = grid['DIMENS']
280
281    # extract only the columns of the face we want to show. the shape was
282    # originally (nelem, 8). copy is necessary to get a contiguous array
283    # for VTK
284    log.debug("Generating cells")
285    src['cells'] = elem_vtcs_ndcs(nk, nj, ni)[:, face].copy()
286    return src
287
288
289def cell_filter(grid, func):
290    """Create a filter for a specific grid.
291
292    :param grid:   Grid structure that should be filtered
293    :type  grid:   dict
294    :param func:   Lambda function that takes i, j, k indices (one-based)
295                   and returns boolean for whether the cells should be
296                   included or not. The function should be vectorized.
297    :type  func:   Callable[(int, int, int), bool]
298
299    >>> layr = cell_filter (grid, lamba i, j, k: np.greater_equal (k, 56))
300    """
301    # extract the dimensions; notice that the i-dimensions is first
302    # in this list, since it is directly from the grid
303    ni, nj, nk = grid['DIMENS']  # pylint: disable=invalid-name
304
305    # setup a grid of the i, j and k address for each of the cells;
306    # notice that now the order of the dimensions are swapped so that
307    # they fit with the memory layout of a loaded Numpy array
308    kji = np.mgrid[1:(nk+1), 1:(nj+1), 1:(ni+1)]
309
310    # call the filter function on all these addresses, and simply
311    # return the boolean array of those
312    filter_flags = func(kji[2], kji[1], kji[0])
313    masked = filter_flags.astype(np.bool)
314
315    # get the mask of active cells, and combine this with the masked
316    # cells from the filter, giving us a flag for all visible nodes
317    active = grid['ACTNUM']
318    visible = np.logical_and(active, masked)
319
320    return visible
321
322
323def face_coords(grid):
324    """Get (x, y, z) coordinates for each corner-point.
325
326    :param grid:  Cornerpoint-grid
327    :type grid:   dict
328    :returns:     Coordinate values for each corner. Use the Face enum to index
329                  the first dimension, k, j, i coordinates to index the next
330                  three, and Dim enum to index the last dimension.
331                  Note that the first point in a face is not necessarily the
332                  point that is closest to the origin of the grid.
333    :rtype:       :class:`numpy.ndarray`, shape = (8, nk, nj, ni, 3)
334
335    >>> import numpy as np
336    >>> import pyresito.io.ecl as ecl
337    >>> import pyresito.grid.cornerpoint as cp
338    >>> case = ecl.EclipseCase ("FOO")
339    >>> coord_fijkd = cp.face_coords (case.grid ())
340    >>> # get the midpoint of the upper face in each cell
341    >>> up_mid = np.average (coord_fijkd [cp.Face.UP, :, :, :, :], axis = 0)
342    """
343    # get coordinates in native corner-point format
344    xyz = corner_coordinates(grid['COORD'], grid['ZCORN'])
345
346    # get extent of grid
347    ni, nj, nk = grid['DIMENS']
348
349    # break up the latter dimension into a many-dimensional hypercube; this
350    # is an inexpensive operation since it won't have to reshuffle the memory
351    xyz = np.reshape(xyz, (3, nk, 2, nj, 2, ni, 2))
352
353    # move all the dimensions local within an element to the front
354    xyz = np.transpose(xyz, (2, 4, 6, 1, 3, 5, 0))
355
356    # coalesce the local dimensions so that they end up in one dimension
357    # which can be indexed by the Face enum to extract a certain field
358    xyz = np.reshape(xyz, (8, nk, nj, ni, 3))
359
360    return xyz
361
362
363def horizon(grid, layer=0, top=True):
364    """Extract the points that are in a certain horizon and average them so
365    that the result is per cell and not per pillar.
366
367    :param grid:   Grid structure
368    :type grid:    dict
369    :param layer:  The K index of the horizon. Default is the layer at the top.
370    :type layer:   int
371    :param top:    Whether the top face should be exported. If this is False,
372                   then the bottom face is exported instead.
373    :type top:     bool
374    :returns:      Depth at the specified horizon for each cell center
375    :rtype:        :class:`numpy.ma.array`, shape = (nk, nj, ni),
376                   dtype = numpy.float64
377    """
378    # get the dimensions of the grid. notice that this array is always given
379    # in Fortran order (since Eclipse works that way) and not in native order
380    num_i, num_j, _ = grid['DIMENS']
381
382    # zcorn is nk * (up,down) * nj * (front,back) * ni * (left,right). we first
383    # project the particular layer and face that we want to work with, ending
384    # up with four corners for each cell.
385    vertical = 0 if top else 1
386    corners = grid['ZCORN'][layer, vertical, :, :, :, :]
387
388    # move the second and fourth dimension (front/back and left/right) to the
389    # back and collapse them into one dimension, meaning that we get
390    # (nj, ni, 4). then take the average across this last dimension, ending up
391    # with an (nj, ni) grid, which is what we return
392    heights = np.average(np.reshape(np.transpose(
393        corners, axes=(0, 2, 1, 3)), (num_j, num_i, 2*2)), axis=2)
394
395    # if there are any inactive elements, then mask out the elevation at that
396    # spot (it is not really a part of the grid)
397    actnum = grid['ACTNUM'][layer, :, :]
398    return ma.array(data=heights, mask=np.logical_not(actnum))
399
400
401def _reduce_corners(plane, red_op):
402    """
403    :param plane:   Z-coordinates for a specific horizon in the grid
404    :type plane:    :class:`numpy.ndarray`, shape = (nj, 2, ni, 2)
405    :param red_op:  How coordinates are merged (reduction operator): either
406                    numpy.minimum or np.maximum, depending on top or bottom
407    :type red_op:   :class:`numpy.ufunc`
408    :returns:       Extremal value of each point around the corners.
409    :rtype:         :class:`numpy.ndarray`, shape = (nj+1, ni+1)
410
411    >> plane = np.reshape (np.arange (36) + 1, (3, 2, 3, 2))
412    >> _reduce_corners (plane, np.minimum)
413    """
414    # constants
415    NORTH = 1
416    SOUTH = 0
417    EAST = 1
418    WEST = 0
419
420    # get dimensions of the domain
421    nj, _, ni, _ = plane.shape
422
423    # allocate output memory for all pillars
424    out = np.empty((nj + 1, ni + 1), dtype=np.float64)
425
426    # views into the original plane, getting a specific corner in every cell.
427    # these views then have dimensions (nj, ni)
428    nw = plane[:, NORTH, :, WEST]
429    ne = plane[:, NORTH, :, EAST]
430    sw = plane[:, SOUTH, :, WEST]
431    se = plane[:, SOUTH, :, EAST]
432
433    # four corners of the output grid only have a single value. the views have
434    # indices 0..n-1, the output array have indices 0..n
435    out[0,  0] = sw[0,    0]
436    out[nj,  0] = nw[nj-1,    0]
437    out[0, ni] = se[0, ni-1]
438    out[nj, ni] = ne[nj-1, ni-1]
439
440    # four edges contains values from two windows; first two horizontal edges,
441    # (varying along i), then two vertical edges (varying along j). recall that
442    # *ranges* in Python have inclusive start, exclusive end.
443    out[0, 1:ni] = red_op(sw[0, 1:ni], se[0, 0:ni-1])
444    out[nj, 1:ni] = red_op(nw[nj-1, 1:ni], ne[nj-1, 0:ni-1])
445    out[1:nj,    0] = red_op(nw[0:nj-1,    0], sw[1:nj,      0])
446    out[1:nj,   ni] = red_op(ne[0:nj-1, ni-1], se[1:nj,   ni-1])
447
448    # interior have two nested reduction operations
449    out[1:nj, 1:ni] = red_op(
450        red_op(nw[0:nj-1,   1:ni], sw[1:nj,   1:ni]),
451        red_op(ne[0:nj-1, 0:ni-1], se[1:nj, 0:ni-1]))
452
453    return out
454
455
456def horizon_pillars(grid, layer=0, top=True):
457    """Extract heights where a horizon crosses the pillars.
458
459    :param grid:  Grid structure, containing both COORD and ZCORN properties.
460    :type grid:   dict
461    :param layer:  The K index of the horizon. Default is the layer at the top.
462    :type layer:   int
463    :param top:    Whether the top face should be exported. If this false is
464                   False, then the bottom face is exported instead.
465    :type top:     bool
466    :returns:     Heights of the horizon at each pillar, in the same format as
467                  the COORD matrix.
468    :rtype:       :class:`numpy.ndarray`, shape = (nj+1, ni+1, 2, 3)
469    """
470    # NumPy's average operations finds the average within arrays, but we want
471    # an operation that can take the average across two arrays
472    def avg_op(a, b):
473        return 0.5 * (a + b)
474
475    # zcorn is nk * (up,down) * nj * (front,back) * ni * (left,right). we first
476    # project the particular layer and face that we want to work with, ending
477    # up with four corners for each cell. take the average of the heights that
478    # are hinged up on each pillar to get a "smooth" surface. (this will not be
479    # the same surface that the cornerpoint grid was created from).
480    vertical = 0 if top else 1
481    points = _reduce_corners(grid['ZCORN'][layer, vertical, :, :, :, :],
482                             avg_op)
483    return points
484
485
486def snugfit(grid):
487    """Create coordinate pillars that matches exactly the top and bottom
488    horizon of the grid cells.
489
490    This version assumes that the pillars in the grid are all strictly
491    vertical, i.e. the x- and y-coordinates will not be changed.
492
493    :param grid:  Grid structure, containing both COORD and ZCORN properties.
494    :type grid:   dict
495    :returns:     Pillar coordinates, in the same format as the COORD matrix.
496                  Note that a new matrix is returned, the original grid is not
497                  altered/updated.
498    :rtype:       :class:`numpy.ndarray`, shape = (nj+1, ni+1, 2, 3)
499    """
500    # placement in the plane of each pillar is unchanged
501    x = grid['COORD'][:, :, 0, 0]
502    y = grid['COORD'][:, :, 0, 1]
503
504    # get implicit dimension of grid
505    njp1, nip1 = x.shape
506
507    # get new top and bottom from the z-coordinates; smaller z is shallower
508    top = _reduce_corners(grid['ZCORN'][0, 0, :, :, :, :], np.minimum)
509    btm = _reduce_corners(grid['ZCORN'][-1, 1, :, :, :, :], np.maximum)
510
511    # combine coordinates to a pillar grid with new z-coordinates
512    coord = np.reshape(np.dstack((
513        x.ravel(), y.ravel(), top.ravel(),
514        x.ravel(), y.ravel(), btm.ravel())), (njp1, nip1, 2, 3))
515
516    return coord
517
518
519def bounding_box(corn, filtr):
520    """\
521    Bounding box of the grid. This function will always assume that the grid
522    is aligned with the geographical axes.
523
524    :param corn:    Coordinate values for each corner. This matrix can be
525                    constructed with the `corner_coordinates` function.
526    :type  corn:    :class:`numpy.ndarray`, shape=(3, nk*2*nj*2*ni*2)
527    :param filtr:   Active corners; use scatter of ACTNUM if no filtering
528    :type  filtr:   :class:`numpy.ndarray`, shape = (nk, 2, nj, 2, ni, 2),
529                    dtype = numpy.bool
530    :return:        Bottom front right corner and top back left corner. Since
531                    the matrix is returned with C ordering, it is specified the
532                    opposite way of what is natural for mathematical matrices.
533    :rtype:         :class:`numpy.ndarray`, shape = (2, 3), dtype = np.float64
534
535    >>> corn = corner_coordinates(grid['COORD'], grid['ZCORN'])
536    >>> bbox = bounding_box(corn, scatter(grid['ACTNUM']))
537    >>> p, q = bbox[0, :], bbox[1, :]
538    >>> diag = np.sqrt(np.dot(q - p, q - p))
539    """
540    # scatter the filter mask from each cell to its corners
541    mask = np.logical_not(filtr).ravel()
542
543    # build arrays of each separate coordinate based on its inclusion
544    x_vals = ma.array(data=corn[0, :], mask=mask)
545    y_vals = ma.array(data=corn[1, :], mask=mask)
546    z_vals = ma.array(data=corn[2, :], mask=mask)
547
548    # construct bounding box that includes all coordinates in visible grid
549    bbox = np.array([[np.min(x_vals), np.min(y_vals), np.min(z_vals)],
550                     [np.max(x_vals), np.max(y_vals), np.max(z_vals)]],
551                    dtype=corn.dtype)
552    return bbox
553
554
555def mass_center(corn, filtr):
556    """\
557    Mass center of the grid. This function will always assume that the density
558    is equal throughout the field.
559
560    :param corn:    Coordinate values for each corner. This matrix can be
561                    constructed with the `corner_coordinates` function.
562    :type  corn:    :class:`numpy.ndarray`, shape=(3, nk*2*nj*2*ni*2)
563    :param filtr:   Active corners; use scatter of ACTNUM if no filtering
564    :type  filtr:   :class:`numpy.ndarray`, shape = (nk, 2, nj, 2, ni, 2),
565                    dtype = numpy.bool
566    :return:        Center of mass. This should be the focal point of the grid.
567    :rtype:         :class:`numpy.ndarray`, shape = (3,), dtype = np.float64
568
569    >>> corn = cp.corner_coordinates(grid['COORD'], grid['ZCORN'])
570    >>> focal_point = cp.mass_center(corn, cp.scatter(grid['ACTNUM']))
571    """
572    # scatter the filter mask from each cell to its corners
573    mask = np.logical_not(filtr).ravel()
574
575    # build arrays of each separate coordinate based on its inclusion
576    x_vals = ma.array(data=corn[0, :], mask=mask)
577    y_vals = ma.array(data=corn[1, :], mask=mask)
578    z_vals = ma.array(data=corn[2, :], mask=mask)
579
580    # use this as a coarse approximation to find the element center
581    center = np.array([np.average(x_vals),
582                       np.average(y_vals),
583                       np.average(z_vals)], dtype=corn.dtype)
584    return center