Merge pull request #826 from bp/gcs_grid_face_arrays_dtype

dtype exposed in GridConnectionSet.grid_face_arrays()

resqpy/fault/_grid_connection_set.py

@@ -2043,7 +2043,8 @@ def grid_face_arrays(self,
                          feature_index = None,
                          active_only = True,
                          lazy = False,
-                         baffle_uuid = None):
+                         baffle_uuid = None,
+                         dtype = None):
         """Creates a triplet of grid face numpy arrays populated from a property for this gcs.
 
         arguments:
@@ -2061,6 +2062,8 @@ def grid_face_arrays(self,
             baffle_uuid (uuid, optional): if present, the uuid of a discrete (bool) property
                 of the gcs holding baffle flags; where True the output face value is set
                 to zero regardless of the main property value
+            dtype (type or str, optional): the element type for the returned arrays; defaults
+                to float for continuous properties or int for discrete properties; see notes
 
         returns:
             triple numpy arrays: identifying the K, J & I direction grid face property values;
@@ -2070,26 +2073,29 @@ def grid_face_arrays(self,
             can only be used on single grid gcs; gcs property must have indexable of faces;
             at present generates grid properties with indexable 'faces' not 'faces per cell',
             which might not be appropriate for grids with split pillars (structural faults);
-            points properties not currently supported; count must be 1
+            points properties not currently supported; count must be 1;
+            if the property is a boolean array and may have been written to hdf5 using packing,
+            then the dtype argument must be set to bool or np.uint8 to ensure unpacking
         """
 
         assert self.number_of_grids() == 1
         (nk, nj, ni) = self.grid_list[0].extent_kji
         active_mask = None
         if active_only:
             pc = self.extract_property_collection()
-            active_mask = pc.single_array_ref(property_kind = 'active')
+            active_mask = pc.single_array_ref(property_kind = 'active', dtype = bool)
             if active_mask is not None:
                 assert active_mask.shape == (self.count,)
         gcs_prop = rqp.Property(self.model, uuid = property_uuid)
         assert gcs_prop is not None
         assert bu.matching_uuids(gcs_prop.collection.support_uuid, self.uuid)
         assert gcs_prop.count() == 1
         assert not gcs_prop.is_points()
-        dtype = float if gcs_prop.is_continuous() else int
+        if dtype is None:
+            dtype = float if gcs_prop.is_continuous() else int
         if default_value is None:
             default_value = -1 if dtype is int else np.nan
-        gcs_prop_array = gcs_prop.array_ref()
+        gcs_prop_array = gcs_prop.array_ref(dtype = dtype)
         assert gcs_prop_array.shape == (self.count,)
         log.debug(f'preparing grid face arrays from gcs property: {gcs_prop.title}; from gcs:{self.title}')
 

tests/unit_tests/fault/test_fault_connection_set.py

@@ -1715,3 +1715,174 @@ def test_feature_index(tmp_path):
     assert np.count_nonzero(j_faces >= 0) == 12
     assert np.all(i_faces == -1)
     assert np.all(j_faces < 112)
+
+
+def test_feature_index_bool_pack(tmp_path):
+    epc = os.path.join(tmp_path, 'gcs_with_packed_bool.epc')
+    model = rq.new_model(epc)
+    crs = rqc.Crs(model)
+    crs.create_xml()
+    grid = grr.RegularGrid(model, extent_kji = (4, 5, 6), crs_uuid = crs.uuid)
+    grid.create_xml()
+    gcs_0 = rqf.GridConnectionSet(model, find_properties = False, grid = grid, title = 'fault zero gcs')
+    gcs_0.set_pairs_from_kelp(
+        kelp_0 = [(2, 0), (2, 1), (2, 2)],  # J face columns
+        kelp_1 = [],  # I face columns
+        feature_name = 'fault zero',
+        create_organizing_objects_where_needed = True)
+    gcs_0.write_hdf5()
+    gcs_0.create_xml()
+    pc = gcs_0.extract_property_collection()
+    a = np.array([False, True] * 6, dtype = bool)
+    pc.add_cached_array_to_imported_list(a,
+                                         'test',
+                                         'this is a pack boolean',
+                                         discrete = True,
+                                         property_kind = 'flag',
+                                         indexable_element = 'faces')
+    pc.write_hdf5_for_imported_list(dtype = np.uint8, use_pack = True)
+    uuids = pc.create_xml_for_imported_list_and_add_parts_to_model()
+    assert len(uuids) == 1
+    kiwi_uuids = uuids
+    gcs_1 = rqf.GridConnectionSet(model, find_properties = False, grid = grid, title = 'fault one gcs')
+    gcs_1.set_pairs_from_kelp(
+        kelp_0 = [],  # J face columns
+        kelp_1 = [(2, 2)],  # I face columns
+        feature_name = 'fault one',
+        create_organizing_objects_where_needed = True)
+    gcs_1.write_hdf5()
+    gcs_1.create_xml()
+    pc = gcs_1.extract_property_collection()
+    a = np.array((True, False, True, True), dtype = bool)
+    pc.add_cached_array_to_imported_list(a,
+                                         'test',
+                                         'this is a single element packed bool!',
+                                         discrete = True,
+                                         property_kind = 'flag',
+                                         indexable_element = 'faces')
+    pc.write_hdf5_for_imported_list(dtype = np.uint8, use_pack = True)
+    uuids = pc.create_xml_for_imported_list_and_add_parts_to_model()
+    assert len(uuids) == 1
+    kiwi_uuids += uuids
+    gcs_2 = rqf.GridConnectionSet(model, find_properties = False, grid = grid, title = 'fault two gcs')
+    gcs_2.set_pairs_from_kelp(
+        kelp_0 = [(1, 3), (1, 4), (1, 5)],  # J face columns
+        kelp_1 = [],  # I face columns
+        feature_name = 'fault two',
+        create_organizing_objects_where_needed = True)
+    gcs_2.write_hdf5()
+    gcs_2.create_xml()
+    pc = gcs_2.extract_property_collection()
+    a = np.array([True] * 8 + [False] * 4, dtype = bool)
+    pc.add_cached_array_to_imported_list(a,
+                                         'test',
+                                         'this is another packed bool',
+                                         discrete = True,
+                                         property_kind = 'flag',
+                                         indexable_element = 'faces')
+    pc.write_hdf5_for_imported_list(dtype = np.uint8, use_pack = True)
+    uuids = pc.create_xml_for_imported_list_and_add_parts_to_model()
+    assert len(uuids) == 1
+    kiwi_uuids += uuids
+    model.store_epc()
+    # check that fault interpretations have been generated
+    fi_titles = model.titles(obj_type = 'FaultInterpretation')
+    assert len(fi_titles) == 3
+    assert set(fi_titles) == set(['fault zero', 'fault one', 'fault two'])
+    # check one of the single feature grid connexion sets
+    fi_one_uuid = model.uuid(obj_type = 'FaultInterpretation', title = 'fault one')
+    assert fi_one_uuid is not None
+    gcs_one_uuid = model.uuid(obj_type = 'GridConnectionSetRepresentation', related_uuid = fi_one_uuid)
+    assert gcs_one_uuid is not None
+    assert bu.matching_uuids(gcs_one_uuid, gcs_1.uuid)
+    gcs_reload = rqf.GridConnectionSet(model, uuid = gcs_one_uuid)
+    assert gcs_reload is not None
+    assert gcs_reload.number_of_grids() == 1
+    assert gcs_reload.number_of_features() == 1
+    assert gcs_reload.title == 'fault one gcs'
+    assert gcs_reload.list_of_feature_names(strip = False) == ['fault one']
+    # make a composite grid connexion set with 3 features
+    gcs = rqf.GridConnectionSet.from_gcs_uuid_list(model,
+                                                   source_model = model,
+                                                   gcs_uuid_list = [gcs_0.uuid, gcs_1.uuid, gcs_2.uuid],
+                                                   title = 'multi feature fault system',
+                                                   gcs_property_uuid_list_of_lists = [kiwi_uuids])
+    model.store_epc()
+    # above class method includes writing hdf5 and creating xml
+    # reopen the composite grid connexion set
+    gcs = rqf.GridConnectionSet(model, uuid = gcs.uuid)
+    assert gcs is not None
+    gcs.cache_arrays()
+    assert gcs.number_of_grids() == 1
+    assert gcs.number_of_features() == 3
+    assert gcs.count == 28
+    # check that a property exists on the composite gcs
+    pc = gcs.extract_property_collection()
+    # a = np.concatenate((np.arange(12, dtype = int) + 100, np.arange(4, dtype = int) + 200, np.arange(12, dtype = int) + 300))
+    # pc.add_cached_array_to_imported_list(a, 'test', 'this is the property of a gcs', discrete = True, null_value = -1,
+    #                                      property_kind = 'kiwi', indexable_element = 'faces')
+    # pc.write_hdf5_for_imported_list()
+    # uuids = pc.create_xml_for_imported_list_and_add_parts_to_model()
+    # assert len(uuids) == 1
+    assert pc.number_of_parts() == 1
+    kiwi_uuid = model.uuid(obj_type = 'DiscreteProperty', related_uuid = gcs.uuid)
+    assert kiwi_uuid is not None
+    # check features
+    assert gcs.list_of_feature_names(strip = False) == ['fault zero', 'fault one', 'fault two']
+    assert gcs.list_of_fault_names(strip = False) == ['fault zero', 'fault one', 'fault two']
+    index, fi_uuid = gcs.feature_index_and_uuid_for_fault_name('fault two')
+    assert index == 2
+    assert fi_uuid is not None
+    fi = rqo.FaultInterpretation(model, uuid = fi_uuid)
+    assert fi is not None
+    assert fi.title == 'fault two'
+    assert gcs.feature_name_for_feature_index(0, strip = False) == 'fault zero'
+    assert gcs.fault_name_for_feature_index(1, strip = False) == 'fault one'
+    # check a few faces are associated with expected features
+    assert gcs.feature_index_for_cell_face((1, 3, 4), 1, 0) is None
+    assert gcs.feature_index_for_cell_face((2, 2, 2), 2, 1) == 1
+    assert gcs.feature_index_for_cell_face((2, 2, 2), 1, 1) == 0
+    assert gcs.feature_index_for_cell_face((0, 2, 4), 1, 0) == 2
+    # check faces associated with each feature
+    k_faces, j_faces, i_faces = gcs.grid_face_arrays(kiwi_uuid,
+                                                     default_value = False,
+                                                     feature_index = 2,
+                                                     active_only = False,
+                                                     lazy = False,
+                                                     baffle_uuid = None,
+                                                     dtype = bool)
+    assert k_faces is not None and j_faces is not None and i_faces is not None
+    assert k_faces.shape == (5, 5, 6)
+    assert j_faces.shape == (4, 6, 6)
+    assert i_faces.shape == (4, 5, 7)
+    assert not np.any(k_faces)
+    assert np.count_nonzero(j_faces) == 8
+    assert not np.any(i_faces)
+    k_faces, j_faces, i_faces = gcs.grid_face_arrays(kiwi_uuid,
+                                                     default_value = False,
+                                                     feature_index = 1,
+                                                     active_only = False,
+                                                     lazy = False,
+                                                     baffle_uuid = None,
+                                                     dtype = bool)
+    assert k_faces is not None and j_faces is not None and i_faces is not None
+    assert k_faces.shape == (5, 5, 6)
+    assert j_faces.shape == (4, 6, 6)
+    assert i_faces.shape == (4, 5, 7)
+    assert not np.any(k_faces)
+    assert not np.any(j_faces)
+    assert np.count_nonzero(i_faces) == 3
+    k_faces, j_faces, i_faces = gcs.grid_face_arrays(kiwi_uuid,
+                                                     default_value = False,
+                                                     feature_index = 0,
+                                                     active_only = False,
+                                                     lazy = False,
+                                                     baffle_uuid = None,
+                                                     dtype = bool)
+    assert k_faces is not None and j_faces is not None and i_faces is not None
+    assert k_faces.shape == (5, 5, 6)
+    assert j_faces.shape == (4, 6, 6)
+    assert i_faces.shape == (4, 5, 7)
+    assert not np.any(k_faces)
+    assert np.count_nonzero(j_faces) == 6
+    assert not np.any(i_faces)