Update CO2 data generation script (#144)

* Update CO2 data generation script

Two main tasks:
- Update to reflect the new format for the CO2 containement table
- Fix issue with reproducibility. Seed should now be set correctly so
  that repeated runs produces the same result.

The latter fix ensures that future changes should not lead to changes in
more files than necessary.

* Update CO2 data generation script, 2

* CO2 data: Add well picks file.

---------

Co-authored-by: Vegard Berg Kvernelv <[email protected]>

01_drogon_ccs/generate-test-data.py

@@ -3,8 +3,9 @@
 generation is done in an ad-hoc and pragmatic manner. The primary purpose is to adhere to
 the folder structure, file formats and naming conventions expected by the plugin.
 
-The script only depends on `input/leakage_boundary.csv` and fault polygon data from
-01_drogon_ahm, all of which are copied to the underlying realization folders.
+The script only depends on `input/leakage_boundary.csv`, `input/hazardous_boundary.csv`
+and fault polygon data from 01_drogon_ahm, all of which are copied to the underlying
+realization folders.
 """
 import datetime
 import pathlib
@@ -45,7 +46,7 @@ def __hash__(self) -> int:
         # NBNB! Ignores x_lin and y_lin. Should not matter as long as the script is set
         # up to always use the same x_lin and y_lin for everything. Pragmatic, but not
         # elegant...
-        return hash(''.join(p.wkt for p in self.polys))
+        return hash("".join(p.wkt for p in self.polys))
 
     def __eq__(self, other) -> bool:
         return hash(self) == hash(other)
@@ -55,11 +56,11 @@ def __eq__(self, other) -> bool:
 def _map_polygons(mp):
     # Determine poly-containment
     res = resolution(mp.x_lin, mp.y_lin)
-    xx, yy = np.meshgrid(mp.x_lin, mp.y_lin, indexing='ij')
+    xx, yy = np.meshgrid(mp.x_lin, mp.y_lin, indexing="ij")
     pts = pygeos.points(xx.flatten(), yy.flatten())
     tree = pygeos.STRtree(pts)
     polys = pygeos.multipolygons([pygeos.from_shapely(p) for p in mp.polys])
-    close = tree.query(polys, 'dwithin', res)
+    close = tree.query(polys, "dwithin", res)
     ij = np.unravel_index(close, shape=xx.shape)
     contained = np.zeros_like(xx, dtype=bool)
     contained[ij] = 1
@@ -71,12 +72,7 @@ def map_polygons(x_lin, y_lin, polys: List[sg.Polygon]):
 
 
 def cost_map(
-    polygon_map: np.ndarray,
-    res: float,
-    distortion_range,
-    min_cost,
-    poly_penalty,
-    seed
+    polygon_map: np.ndarray, res: float, distortion_range, min_cost, poly_penalty, seed
 ):
     gen = np.random.RandomState(seed)
     field = gen.normal(size=polygon_map.shape)
@@ -85,8 +81,7 @@ def cost_map(
     dist_field = np.exp(dist_field)
     dist_field = dist_field / np.std(dist_field) + min_cost
     poly_field = scipy.ndimage.gaussian_filter(
-        1 + polygon_map * float(poly_penalty),
-        sigma=distortion_range / 10
+        1 + polygon_map * float(poly_penalty), sigma=distortion_range / 10
     )
     mult = dist_field * poly_field
     return mult * res
@@ -142,7 +137,8 @@ def simulate_migration_over_time(
     **emulate_kwargs,
 ):
     results = {
-        t: scaling * emulate_migration(
+        t: scaling
+        * emulate_migration(
             **emulate_kwargs,
             migration_dist=init_mig_dist + i * init_mig_dist / 15,
         )
@@ -164,9 +160,9 @@ def generate_maps(output_dir, surface_name, time_steps, init_mig_dist, **kwargs)
     if not output_dir.is_dir():
         output_dir.mkdir(parents=True)
     sgas = simulate_migration_over_time(time_steps, 0.6, init_mig_dist, **kwargs)
-    kwargs['seed'] += 999
+    kwargs["seed"] += 999
     amfg = simulate_migration_over_time(time_steps, 1e-4, init_mig_dist, **kwargs)
-    template = kwargs['template']
+    template = kwargs["template"]
     # SGAS
     for t, s in sgas.items():
         surf = array_to_xtgeo(template, s, 1e-4)
@@ -176,10 +172,7 @@ def generate_maps(output_dir, surface_name, time_steps, init_mig_dist, **kwargs)
         surf = array_to_xtgeo(template, s, 1e-8)
         surf.to_file(output_dir / f"{surface_name}--max_AMFG--{t}.gri")
     # Migration Time
-    mtime_all = [
-        np.where(s > 1e-2, float(t[:4]), np.inf)
-        for t, s in sgas.items()
-    ]
+    mtime_all = [np.where(s > 1e-2, float(t[:4]), np.inf) for t, s in sgas.items()]
     mtime = np.min(mtime_all, axis=0)
     mtime -= np.min(mtime)
     surf = array_to_xtgeo(template, mtime, -np.inf, 999999999)
@@ -188,19 +181,52 @@ def generate_maps(output_dir, surface_name, time_steps, init_mig_dist, **kwargs)
 
 
 def simulate_containment(
-    tmpl: xtgeo.RegularSurface, saturation, polygon: sg.Polygon, seed
+    tmpl: xtgeo.RegularSurface,
+    saturation,
+    containment_boundary: sg.Polygon,
+    hazardous_boundary: sg.Polygon,
+    seed,
+    calc_volume: bool = False,
 ):
     x_lin = tmpl.xmin + np.arange(tmpl.ncol) * tmpl.xinc
     y_lin = tmpl.ymin + np.arange(tmpl.nrow) * tmpl.yinc
-    poly_map = map_polygons(x_lin, y_lin, [polygon])
+    poly_map_con = map_polygons(x_lin, y_lin, [containment_boundary])
+    poly_map_haz = map_polygons(x_lin, y_lin, [hazardous_boundary])
     gen = np.random.RandomState(seed)
-    volumes = gen.uniform(7, 13, size=poly_map.shape) * resolution(x_lin, y_lin) ** 2
-    poro = 0.3
-    density = 700
-    mass = poro * volumes * density * saturation
-    within = mass[poly_map].sum()
-    outside = mass[~poly_map].sum()
-    return within, outside
+    volumes = (
+        gen.uniform(7, 13, size=poly_map_con.shape) * resolution(x_lin, y_lin) ** 2
+    )
+    if calc_volume:
+        prop = volumes * saturation
+    else:  # Calculate CO2 mass
+        poro = 0.3
+        density = 700
+        prop = poro * volumes * density * saturation
+    is_contained, is_outside = [], []
+    for a, b in zip(poly_map_con, poly_map_haz):
+        is_contained.append([x if not y else False for x, y in zip(a, b)])
+        is_outside.append([not x and not y for x, y in zip(a, b)])
+    is_contained = np.array(is_contained)
+    is_outside = np.array(is_outside)
+    contained = prop[is_contained].sum()
+    outside = prop[is_outside].sum()
+    hazardous = prop[poly_map_haz].sum()
+
+    # Simulate fractions
+    aqu_contained = contained * gen.uniform(high=0.1)
+    gas_contained = contained - aqu_contained
+    aqu_outside = outside * gen.uniform(high=0.1)
+    gas_outside = outside - aqu_outside
+    aqu_hazardous = hazardous * gen.uniform(high=0.1)
+    gas_hazardous = hazardous - aqu_hazardous
+    return (
+        aqu_contained,
+        gas_contained,
+        aqu_outside,
+        gas_outside,
+        aqu_hazardous,
+        gas_hazardous,
+    )
 
 
 def setup_ensemble_folders(ens_root, input_folder, polygons_folder):
@@ -215,9 +241,9 @@ def setup_ensemble_folders(ens_root, input_folder, polygons_folder):
     for f in polygons_folder.glob("*gl_faultlines_extract_postprocess.pol"):
         shutil.copy(f, res_root / "polygons")
     shutil.copy(input_folder / "leakage_boundary.csv", res_root / "polygons")
+    shutil.copy(input_folder / "hazardous_boundary.csv", res_root / "polygons")
     # Write dummy OK and STATUS files
-    t = datetime.datetime.now()
-    t = datetime.datetime.strftime(t, "%H:%M:%S")
+    t = "13:42:37"
     with open(ens_root / "iter-0" / "OK", "w") as f:
         f.write(f"All jobs complete {t}\n")
     with open(ens_root / "iter-0" / "STATUS", "w") as f:
@@ -226,14 +252,53 @@ def setup_ensemble_folders(ens_root, input_folder, polygons_folder):
     return res_root
 
 
-def main(ens_root, input_folder, polygons_folder):
+def generate_date_table_entry(
+    surface_template: xtgeo.RegularSurface,
+    saturation: xtgeo.RegularSurface,
+    containment_boundary: sg.Polygon,
+    hazardous_boundary: sg.Polygon,
+    seed: int,
+    calc_volume: bool = False,
+):
+    ca, cg, oa, og, ha, hg = simulate_containment(
+        surface_template,
+        saturation,
+        containment_boundary,
+        hazardous_boundary,
+        seed,
+        calc_volume,
+    )
+    return {
+        "total": ca + cg + oa + og + ha + hg,
+        "total_contained": ca + cg,
+        "total_outside": oa + og,
+        "total_hazardous": ha + hg,
+        "total_gas": cg + og + hg,
+        "total_aqueous": ca + oa + ha,
+        "gas_contained": cg,
+        "aqueous_contained": ca,
+        "gas_outside": og,
+        "aqueous_outside": oa,
+        "gas_hazardous": hg,
+        "aqueous_hazardous": ha,
+    }
+
+
+def main(ens_root, input_folder, polygons_folder, base_seed):
     res_root = setup_ensemble_folders(ens_root, input_folder, polygons_folder)
     polys = read_polylines(
         res_root / "polygons" / "topvolantis--gl_faultlines_extract_postprocess.csv"
     )
-    boundary = sg.Polygon(np.genfromtxt(
-        input_folder / "leakage_boundary.csv", skip_header=1, delimiter=','
-    ))
+    containment_boundary = sg.Polygon(
+        np.genfromtxt(
+            input_folder / "leakage_boundary.csv", skip_header=1, delimiter=","
+        )
+    )
+    hazardous_boundary = sg.Polygon(
+        np.genfromtxt(
+            input_folder / "hazardous_boundary.csv", skip_header=1, delimiter=","
+        )
+    )
     tmpl = xtgeo.RegularSurface(
         ncol=279, nrow=341, xinc=25.0, yinc=25.0, xori=460063.6875, yori=5929551.0
     )
@@ -244,8 +309,8 @@ def main(ens_root, input_folder, polygons_folder):
         "topvolon": 70,
         "toptherys": 70,
     }
-    containments = []
-    base_seed = hash(ens_root) % 2**32
+    mass_containments = []
+    volume_containments = []
     for sn, imd in mig_dists.items():
         sgas, amfg = generate_maps(
             res_root / "maps",
@@ -261,23 +326,53 @@ def main(ens_root, input_folder, polygons_folder):
             cutoff=8,
             seed=base_seed,
         )
-        containments += [
-            [sn, t, *simulate_containment(tmpl, s, boundary, base_seed - 1)]
+        mass_containments += [
+            dict(
+                **generate_date_table_entry(
+                    tmpl,
+                    s,
+                    containment_boundary,
+                    hazardous_boundary,
+                    (base_seed + 1) % 2**32,
+                ),
+                date=f"{t[:4]}-{t[4:6]}-{t[6:8]}",
+            )
+            for t, s in sgas.items()
+        ]
+        volume_containments += [
+            dict(
+                **generate_date_table_entry(
+                    tmpl,
+                    s,
+                    containment_boundary,
+                    hazardous_boundary,
+                    (base_seed + 1) % 2**32,
+                    True,
+                ),
+                date=f"{t[:4]}-{t[4:6]}-{t[6:8]}",
+            )
             for t, s in sgas.items()
         ]
-    df = pd.DataFrame.from_records(
-        containments,
-        columns=['surface', 'date', 'co2_inside', 'co2_outside'],
+    df_mass = pd.DataFrame.from_records(mass_containments)
+    df_plume_volume_actual = pd.DataFrame.from_records(volume_containments)
+    df_mass = df_mass.groupby("date").sum()
+    df_plume_volume_actual = df_plume_volume_actual.groupby("date").sum()
+    df_plume_volume_actual_simple = df_plume_volume_actual * 0.8
+
+    df_mass.to_csv(res_root / "tables/co2_volumes.csv")
+    df_plume_volume_actual.to_csv(res_root / "tables/plume_volume_actual.csv")
+    df_plume_volume_actual_simple.to_csv(
+        res_root / "tables/plume_volume_actual_simple.csv"
     )
-    df.groupby('date').sum().to_csv(res_root / "tables/co2_volumes.csv")
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     for i in tqdm(range(5)):
         current = pathlib.Path(__file__).parent
         main(
             current / f"realization-{i}",
             current / "input",
-            current / "../01_drogon_ahm/realization-0/iter-0/share/results/polygons"
+            current / "../01_drogon_ahm/realization-0/iter-0/share/results/polygons",
+            i + 42,
         )
     print(f"Cache efficiency: {_map_polygons.cache_info()}")

01_drogon_ccs/input/hazardous_boundary.csv

@@ -0,0 +1,7 @@
+x,y,z
+462154,5933414,0
+462381,5932578,0
+462919,5932320,0
+460000,5932320,0
+460000,5933414,0
+462154,5933414,0

01_drogon_ccs/realization-0/iter-0/OK

@@ -1 +1 @@
-All jobs complete 15:01:28
+All jobs complete 13:42:37

01_drogon_ccs/realization-0/iter-0/STATUS

@@ -1,2 +1,2 @@
 Current host                    : fake-host
-MAKE_DIRECTORY                  : 15:01:28 .... 15:01:28
+MAKE_DIRECTORY                  : 13:42:37 .... 13:42:37

01_drogon_ccs/realization-0/iter-0/share/results/polygons/hazardous_boundary.csv

@@ -0,0 +1,7 @@
+x,y,z
+462154,5933414,0
+462381,5932578,0
+462919,5932320,0
+460000,5932320,0
+460000,5933414,0
+462154,5933414,0

01_drogon_ccs/realization-0/iter-0/share/results/tables/co2_volumes.csv

@@ -1,11 +1,11 @@
-date,co2_inside,co2_outside
-20200101,2899574400.826876,203393917.238401
-20400101,3088392407.3065677,283295051.0170409
-20600101,3266875342.2931423,373310702.3053861
-20800101,3436658691.3793783,473462734.4376221
-21000101,3598976832.9913793,585443645.9197086
-21200101,3754781699.6663804,712232902.9234608
-21400101,3904946618.699225,857545511.2494102
-21600101,4050401631.505111,1025315075.6815368
-21800101,4192159976.0840845,1219303266.212956
-22000101,4331281612.097786,1442848607.331028
+date,total,total_contained,total_outside,total_hazardous,total_gas,total_aqueous,gas_contained,aqueous_contained,gas_outside,aqueous_outside,gas_hazardous,aqueous_hazardous
+2020-01-01,1883742663.4364204,1815824246.75565,67827159.99937923,91256.68139133565,1729507479.4153974,154235184.02102315,1661834525.0888603,153989721.66678983,67588134.51346314,239025.48591609587,84819.81307410785,6436.86831722779
+2040-01-01,2128717771.4416158,2007508582.176139,120951697.65309899,257491.6123779341,1958027991.2769053,170689780.16471064,1837263202.7759404,170245379.40019864,120525459.28038555,426238.37271342764,239329.22057935628,18162.391798577813
+2060-01-01,2391343980.7116237,2194584868.498756,196112681.49108526,646430.7217823985,2204497020.922794,186846959.78882974,2008474613.7876635,186110254.71109265,195421572.9589258,691108.5321594508,600834.1762047629,45596.545577635676
+2080-01-01,2668669491.6071606,2373712354.920289,293486495.4529434,1470641.2339275996,2466230463.5840726,202439028.02308798,2172411317.38635,201301037.53393894,292452237.90499216,1034257.5479512026,1366908.2927297605,103732.941197839
+2100-01-01,2957081455.1234307,2542850289.997418,411155502.6047281,3075662.5212844624,2739770919.682467,217310535.44096377,2327205626.6459703,215644663.35144764,409706575.00997,1448927.5947580764,2858718.026526399,216944.49475806364
+2120-01-01,3253247375.0617323,2701222532.4708047,546044885.7701776,5979956.820750314,3021825978.736671,231421396.3250621,2472147220.4308643,229075312.03994045,544120602.8700176,1924282.9001599108,5558155.435788571,421801.3849617432
+2140-01-01,3554720118.1839566,2849021491.033886,694793701.4838388,10904925.666231468,3309893149.838347,244826968.3456092,2607412190.348054,241609300.6858322,692345222.1120039,2448479.3718347833,10135737.378289254,769188.2879422131
+2160-01-01,3860166743.4822946,2987086912.461138,854297811.4806368,18782019.540520243,3602513516.276303,257653227.20599213,2733769069.0265365,253317843.43460163,851287233.571904,3010577.908732805,17457213.677862536,1324805.8626577077
+2180-01-01,4169323360.6266875,3116646424.531534,1021953364.1085075,30723571.986645605,3899249796.265426,270073564.36126137,2852341376.1189995,264305048.4125348,1018351961.669693,3601402.4388144882,28556458.47673352,2167113.509912086
+2200-01-01,4482802529.143365,3239127375.452816,1195726551.5316484,47948602.15890018,4200514686.7066507,282287842.43671453,2964435414.553817,274691960.8989997,1191512765.785566,4213785.746082332,44566506.367267676,3382095.791632504

01_drogon_ccs/realization-0/iter-0/share/results/tables/plume_volume_actual.csv

@@ -0,0 +1,11 @@
+date,total,total_contained,total_outside,total_hazardous,total_gas,total_aqueous,gas_contained,aqueous_contained,gas_outside,aqueous_outside,gas_hazardous,aqueous_hazardous
+2020-01-01,8970203.159221051,8646782.127407858,322986.47618752025,434.5556256730269,8235749.901978082,734453.2572429674,7913497.7385183815,733284.3888894753,321848.2595879198,1138.216599600457,403.90387178146597,30.651753891560904
+2040-01-01,10136751.292579126,9559564.677029233,575960.465014757,1226.1505351330195,9323942.815604312,812808.4769748125,8748872.394171145,810692.2828580886,573930.7584780264,2029.7065367306077,1139.6629551397918,86.48757999322767
+2060-01-01,11387352.289102972,10450404.135708366,933869.9118623108,3078.241532297135,10497604.861537118,889747.4275658561,9564164.827560306,886239.3081480605,930578.9188520277,3290.993010283099,2861.115124784584,217.12640751255077
+2080-01-01,12707949.9600341,11303392.166287092,1397554.7402521118,7003.053494893332,11743954.588495584,963995.3715385144,10344815.797077859,958576.3692092332,1392629.704309487,4925.035942624775,6509.087108236956,493.9663866563762
+2100-01-01,14081340.262492528,12108810.90474961,1957883.345736801,14646.012006116485,13046528.188964128,1034812.073528399,11081931.555457002,1026879.3492926079,1950983.6905236673,6899.655213133698,13612.94298345904,1033.0690226574457
+2120-01-01,15491654.166960632,12862964.440337166,2600213.7417627517,28475.984860715773,14389647.51779367,1102006.6491669624,11772129.621099355,1090834.8192378115,2591050.4898572285,9163.25190552339,26467.406837088423,2008.578023627348
+2140-01-01,16927238.658018842,13566769.00492327,3308541.435637329,51928.21745824508,15761395.951611178,1165842.706407663,12416248.525466925,1150520.479456344,3296882.0100571634,11659.425580165642,48265.41608709169,3662.8013711533954
+2160-01-01,18381746.39753474,14224223.392672084,4068084.8165744613,89438.18828819165,17154826.267982394,1226920.1295523434,13017947.947745409,1206275.4449266745,4053748.7312947814,14336.085279680023,83129.58894220256,6308.599345989086
+2180-01-01,19853920.76488899,14841173.450150168,4866444.590992892,146302.72374593146,18567856.172692504,1286064.5921964832,13582577.981519049,1258595.4686311183,4849295.055569966,17149.535422926132,135983.13560349296,10319.588142438506
+2200-01-01,21346678.710206494,15424416.073584838,5693935.959674514,228326.67694714371,20002450.889079284,1344227.821127212,14116359.116922935,1308056.9566619033,5673870.313264598,20065.64640991586,212221.45889175084,16105.218055392877

01_drogon_ccs/realization-0/iter-0/share/results/tables/plume_volume_actual_simple.csv

@@ -0,0 +1,11 @@
+date,total,total_contained,total_outside,total_hazardous,total_gas,total_aqueous,gas_contained,aqueous_contained,gas_outside,aqueous_outside,gas_hazardous,aqueous_hazardous
+2020-01-01,7176162.527376842,6917425.701926287,258389.18095001622,347.64450053842154,6588599.921582466,587562.6057943739,6330798.190814706,586627.5111115802,257478.60767033586,910.5732796803655,323.1230974251728,24.521403113248724
+2040-01-01,8109401.034063301,7647651.741623387,460768.37201180565,980.9204281064157,7459154.25248345,650246.78157985,6999097.915336916,648553.826286471,459144.6067824211,1623.7652293844862,911.7303641118335,69.19006399458213
+2060-01-01,9109881.831282377,8360323.308566693,747095.9294898487,2462.5932258377084,8398083.889229694,711797.942052685,7651331.862048245,708991.4465184484,744463.1350816223,2632.7944082264794,2288.8920998276676,173.70112601004064
+2080-01-01,10166359.96802728,9042713.733029675,1118043.7922016894,5602.442795914666,9395163.670796467,771196.2972308116,8275852.637662288,766861.0953673866,1114103.7634475897,3940.02875409982,5207.269686589565,395.173109325101
+2100-01-01,11265072.209994024,9687048.723799689,1566306.6765894408,11716.80960489319,10437222.551171303,827849.6588227192,8865545.2443656,821503.4794340864,1560786.9524189339,5519.724170506959,10890.354386767232,826.4552181259565
+2120-01-01,12393323.333568506,10290371.552269734,2080170.9934102015,22780.78788857262,11511718.014234938,881605.31933357,9417703.696879484,872667.8553902493,2072840.3918857828,7330.601524418712,21173.92546967074,1606.8624189018785
+2140-01-01,13541790.926415075,10853415.203938616,2646833.1485098633,41542.57396659607,12609116.761288943,932674.1651261305,9932998.82037354,920416.3835650752,2637505.6080457307,9327.540464132513,38612.33286967335,2930.2410969227167
+2160-01-01,14705397.118027791,11379378.714137668,3254467.853259569,71550.55063055332,13723861.014385916,981536.1036418747,10414358.358196327,965020.3559413396,3242998.9850358255,11468.868223744019,66503.67115376206,5046.87947679127
+2180-01-01,15883136.611911193,11872938.760120135,3893155.672794314,117042.17899674518,14854284.938154005,1028851.6737571866,10866062.38521524,1006876.3749048947,3879436.0444559734,13719.628338340906,108786.50848279438,8255.670513950805
+2200-01-01,17077342.968165196,12339532.85886787,4555148.767739612,182661.34155771497,16001960.711263428,1075382.2569017697,11293087.293538349,1046445.5653295227,4539096.250611679,16052.517127932688,169777.16711340068,12884.174444314303