diff --git a/autotest/regression/test_model_netcdf.py b/autotest/regression/test_model_netcdf.py index 481eb6e05..a0414d506 100644 --- a/autotest/regression/test_model_netcdf.py +++ b/autotest/regression/test_model_netcdf.py @@ -72,26 +72,18 @@ def compare_netcdf_data(base, gen): # coordinates for coordname, da in xrb.coords.items(): - assert coordname in xrg.coords - print(f"NetCDF file check data equivalence for coordinate: {coordname}") - assert np.allclose(xrb.coords[coordname].data, xrg.coords[coordname].data) + compare_netcdf_var(coordname, xrb.coords, xrg.data_vars, xrg.coords) # variables for varname, da in xrb.data_vars.items(): if varname == "projection": continue - compare_netcdf_var( - varname, - xrb.data_vars, - xrg.data_vars, - xrg.coords, - ) + compare_netcdf_var(varname, xrb.data_vars, xrg.data_vars, xrg.coords) def compare_netcdf_var(varname, base_d, gen_d, coord_d, projection=False, update=None): # check variable name - # assert varname in xrg.data_vars assert varname in gen_d or varname in coord_d if varname in gen_d: @@ -101,11 +93,16 @@ def compare_netcdf_var(varname, base_d, gen_d, coord_d, projection=False, update # encodings for e in base_d[varname].encoding: + assert e in var_d[varname].encoding if e.lower() == "source": continue - assert e in var_d[varname].encoding if e == "_FillValue": - assert np.allclose(base_d[varname].encoding[e], var_d[varname].encoding[e]) + if np.isnan(base_d[varname].encoding[e]): + assert np.isnan(var_d[varname].encoding[e]) + else: + assert np.allclose( + base_d[varname].encoding[e], var_d[varname].encoding[e] + ) else: assert base_d[varname].encoding[e] == var_d[varname].encoding[e] @@ -140,39 +137,39 @@ def test_load_gwfsto01(function_tmpdir, example_data_path): }, } ws = function_tmpdir / "ws" - for base_folder, test_info in tests.items(): - data_path = os.path.join(data_path_base, base_folder, test_info["base_sim_dir"]) + for dirname, test in tests.items(): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) # copy example data into working directory - base_model_folder = os.path.join(ws, f"{base_folder}_base") - test_model_folder = os.path.join(ws, f"{base_folder}_test") - shutil.copytree(data_path, base_model_folder) + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) # load example - sim = flopy.mf6.MFSimulation.load(sim_ws=base_model_folder) + sim = flopy.mf6.MFSimulation.load(sim_ws=base_path) gwf = sim.get_model("gwf_sto01") # set simulation path and write simulation - sim.set_sim_path(test_model_folder) - sim.write_simulation(netcdf=test_info["netcdf_type"]) + sim.set_sim_path(test_path) + sim.write_simulation(netcdf=test["netcdf_type"]) # compare generated files gen_files = [ f - for f in os.listdir(test_model_folder) - if os.path.isfile(os.path.join(test_model_folder, f)) + for f in os.listdir(test_path) + if os.path.isfile(os.path.join(test_path, f)) ] base_files = [ f - for f in os.listdir(base_model_folder) - if os.path.isfile(os.path.join(base_model_folder, f)) + for f in os.listdir(base_path) + if os.path.isfile(os.path.join(base_path, f)) ] assert len(gen_files) == len(base_files) for f in base_files: - base = os.path.join(base_model_folder, f) - gen = os.path.join(test_model_folder, f) - if f != test_info["netcdf_output_file"]: + base = os.path.join(base_path, f) + gen = os.path.join(test_path, f) + if f != test["netcdf_output_file"]: # "gwf_sto01.dis.ncf": # TODO wkt string missing on write? with open(base, "r") as file1, open(gen, "r") as file2: # Skip first line @@ -222,16 +219,16 @@ def test_update_gwfsto01(function_tmpdir, example_data_path): } ws = function_tmpdir / "ws" - for base_folder, test_info in tests.items(): - data_path = os.path.join(data_path_base, base_folder, test_info["base_sim_dir"]) + for dirname, test in tests.items(): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) # copy example data into working directory - base_model_folder = os.path.join(ws, f"{base_folder}_base") - test_model_folder = os.path.join(ws, f"{base_folder}_test") - shutil.copytree(data_path, base_model_folder) + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) # load example - sim = flopy.mf6.MFSimulation.load(sim_ws=base_model_folder) + sim = flopy.mf6.MFSimulation.load(sim_ws=base_path) # get model instance gwf = sim.get_model("gwf_sto01") @@ -244,26 +241,26 @@ def test_update_gwfsto01(function_tmpdir, example_data_path): gwf.ic.strt.set_data(ic_strt) # set simulation path and write simulation - sim.set_sim_path(test_model_folder) - sim.write_simulation(netcdf=test_info["netcdf_type"]) + sim.set_sim_path(test_path) + sim.write_simulation(netcdf=test["netcdf_type"]) # compare generated files gen_files = [ f - for f in os.listdir(test_model_folder) - if os.path.isfile(os.path.join(test_model_folder, f)) + for f in os.listdir(test_path) + if os.path.isfile(os.path.join(test_path, f)) ] base_files = [ f - for f in os.listdir(base_model_folder) - if os.path.isfile(os.path.join(base_model_folder, f)) + for f in os.listdir(base_path) + if os.path.isfile(os.path.join(base_path, f)) ] assert len(gen_files) == len(base_files) for f in base_files: - base = os.path.join(base_model_folder, f) - gen = os.path.join(test_model_folder, f) - if f != test_info["netcdf_output_file"]: + base = os.path.join(base_path, f) + gen = os.path.join(test_path, f) + if f != test["netcdf_output_file"]: # "gwf_sto01.dis.ncf": # TODO wkt string missing on write? with open(base, "r") as file1, open(gen, "r") as file2: # Skip first line @@ -369,13 +366,13 @@ def test_create_gwfsto01(function_tmpdir, example_data_path): # build ws = function_tmpdir / "ws" - for idx, (base_folder, test_info) in enumerate(tests.items()): - data_path = os.path.join(data_path_base, base_folder, test_info["base_sim_dir"]) + for idx, (dirname, test) in enumerate(tests.items()): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) # copy example data into working directory - base_model_folder = os.path.join(ws, f"{base_folder}_base") - test_model_folder = os.path.join(ws, f"{base_folder}_test") - shutil.copytree(data_path, base_model_folder) + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) # build MODFLOW 6 files sim = flopy.mf6.MFSimulation( @@ -486,26 +483,26 @@ def test_create_gwfsto01(function_tmpdir, example_data_path): ) # set simulation path and write simulation - sim.set_sim_path(test_model_folder) - sim.write_simulation(netcdf=test_info["netcdf_type"]) + sim.set_sim_path(test_path) + sim.write_simulation(netcdf=test["netcdf_type"]) # compare generated files gen_files = [ f - for f in os.listdir(test_model_folder) - if os.path.isfile(os.path.join(test_model_folder, f)) + for f in os.listdir(test_path) + if os.path.isfile(os.path.join(test_path, f)) ] base_files = [ f - for f in os.listdir(base_model_folder) - if os.path.isfile(os.path.join(base_model_folder, f)) + for f in os.listdir(base_path) + if os.path.isfile(os.path.join(base_path, f)) ] assert len(gen_files) == len(base_files) for f in base_files: - base = os.path.join(base_model_folder, f) - gen = os.path.join(test_model_folder, f) - if f != test_info["netcdf_output_file"]: + base = os.path.join(base_path, f) + gen = os.path.join(test_path, f) + if f != test["netcdf_output_file"]: with open(base, "r") as file1, open(gen, "r") as file2: # Skip first line next(file1) @@ -595,14 +592,14 @@ def test_gwfsto01(function_tmpdir, example_data_path): strt_longname = "starting head" ws = function_tmpdir / "ws" - for base_folder, test_info in tests.items(): - data_path = os.path.join(data_path_base, base_folder, test_info["base_sim_dir"]) + for dirname, test in tests.items(): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) # copy example data into working directory - base_model_folder = os.path.join(ws, f"{base_folder}_base") - test_model_folder = os.path.join(ws, f"{base_folder}_test") - shutil.copytree(data_path, base_model_folder) - os.mkdir(test_model_folder) + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) + os.mkdir(test_path) # create discretization dis = flopy.discretization.StructuredGrid( @@ -620,8 +617,8 @@ def test_gwfsto01(function_tmpdir, example_data_path): ds = create_dataset( "gwf6", "gwf_sto01", - test_info["netcdf_type"], - test_info["netcdf_output_file"], + test["netcdf_type"], + test["netcdf_output_file"], dis, ) @@ -651,12 +648,12 @@ def test_gwfsto01(function_tmpdir, example_data_path): ds.create_array("sto", "sy", sy, ["nlay", "nrow", "ncol"], sy_longname) # write to netcdf - ds.write(test_model_folder) + ds.write(test_path) # compare compare_netcdf( - os.path.join(base_model_folder, test_info["netcdf_output_file"]), - os.path.join(test_model_folder, test_info["netcdf_output_file"]), + os.path.join(base_path, test["netcdf_output_file"]), + os.path.join(test_path, test["netcdf_output_file"]), projection=True, ) @@ -671,38 +668,38 @@ def test_load_disv01b(function_tmpdir, example_data_path): }, } ws = function_tmpdir / "ws" - for base_folder, test_info in tests.items(): - data_path = os.path.join(data_path_base, base_folder, test_info["base_sim_dir"]) + for dirname, test in tests.items(): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) # copy example data into working directory - base_model_folder = os.path.join(ws, f"{base_folder}_base") - test_model_folder = os.path.join(ws, f"{base_folder}_test") - shutil.copytree(data_path, base_model_folder) + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) # load example - sim = flopy.mf6.MFSimulation.load(sim_ws=base_model_folder) + sim = flopy.mf6.MFSimulation.load(sim_ws=base_path) # set simulation path and write simulation - sim.set_sim_path(test_model_folder) + sim.set_sim_path(test_path) sim.write_simulation(netcdf="mesh2d") # compare generated files gen_files = [ f - for f in os.listdir(test_model_folder) - if os.path.isfile(os.path.join(test_model_folder, f)) + for f in os.listdir(test_path) + if os.path.isfile(os.path.join(test_path, f)) ] base_files = [ f - for f in os.listdir(base_model_folder) - if os.path.isfile(os.path.join(base_model_folder, f)) + for f in os.listdir(base_path) + if os.path.isfile(os.path.join(base_path, f)) ] assert len(gen_files) == len(base_files) for f in base_files: - base = os.path.join(base_model_folder, f) - gen = os.path.join(test_model_folder, f) - if f != test_info["netcdf_output_file"]: + base = os.path.join(base_path, f) + gen = os.path.join(test_path, f) + if f != test["netcdf_output_file"]: with open(base, "r") as file1, open(gen, "r") as file2: # Skip first line next(file1) @@ -755,16 +752,16 @@ def test_update_disv01b(function_tmpdir, example_data_path): } ws = function_tmpdir / "ws" - for base_folder, test_info in tests.items(): - data_path = os.path.join(data_path_base, base_folder, test_info["base_sim_dir"]) + for dirname, test in tests.items(): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) # copy example data into working directory - base_model_folder = os.path.join(ws, f"{base_folder}_base") - test_model_folder = os.path.join(ws, f"{base_folder}_test") - shutil.copytree(data_path, base_model_folder) + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) # load example - sim = flopy.mf6.MFSimulation.load(sim_ws=base_model_folder) + sim = flopy.mf6.MFSimulation.load(sim_ws=base_path) # get model instance gwf = sim.get_model("disv01b") @@ -777,26 +774,26 @@ def test_update_disv01b(function_tmpdir, example_data_path): gwf.ic.strt.set_data(strt) # set simulation path and write simulation - sim.set_sim_path(test_model_folder) + sim.set_sim_path(test_path) sim.write_simulation(netcdf="mesh2d") # compare generated files gen_files = [ f - for f in os.listdir(test_model_folder) - if os.path.isfile(os.path.join(test_model_folder, f)) + for f in os.listdir(test_path) + if os.path.isfile(os.path.join(test_path, f)) ] base_files = [ f - for f in os.listdir(base_model_folder) - if os.path.isfile(os.path.join(base_model_folder, f)) + for f in os.listdir(base_path) + if os.path.isfile(os.path.join(base_path, f)) ] assert len(gen_files) == len(base_files) for f in base_files: - base = os.path.join(base_model_folder, f) - gen = os.path.join(test_model_folder, f) - if f != test_info["netcdf_output_file"]: + base = os.path.join(base_path, f) + gen = os.path.join(test_path, f) + if f != test["netcdf_output_file"]: with open(base, "r") as file1, open(gen, "r") as file2: # Skip first line next(file1) @@ -836,13 +833,13 @@ def test_create_disv01b(function_tmpdir, example_data_path): # build ws = function_tmpdir / "ws" - for idx, (base_folder, test_info) in enumerate(tests.items()): - data_path = os.path.join(data_path_base, base_folder, test_info["base_sim_dir"]) + for idx, (dirname, test) in enumerate(tests.items()): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) # copy example data into working directory - base_model_folder = os.path.join(ws, f"{base_folder}_base") - test_model_folder = os.path.join(ws, f"{base_folder}_test") - shutil.copytree(data_path, base_model_folder) + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) # create simulation sim = flopy.mf6.MFSimulation( @@ -868,26 +865,26 @@ def test_create_disv01b(function_tmpdir, example_data_path): ) # set path and write simulation - sim.set_sim_path(test_model_folder) - sim.write_simulation(netcdf=test_info["netcdf_type"]) + sim.set_sim_path(test_path) + sim.write_simulation(netcdf=test["netcdf_type"]) # compare generated files gen_files = [ f - for f in os.listdir(test_model_folder) - if os.path.isfile(os.path.join(test_model_folder, f)) + for f in os.listdir(test_path) + if os.path.isfile(os.path.join(test_path, f)) ] base_files = [ f - for f in os.listdir(base_model_folder) - if os.path.isfile(os.path.join(base_model_folder, f)) + for f in os.listdir(base_path) + if os.path.isfile(os.path.join(base_path, f)) ] assert len(gen_files) == len(base_files) for f in base_files: - base = os.path.join(base_model_folder, f) - gen = os.path.join(test_model_folder, f) - if f != test_info["netcdf_output_file"]: + base = os.path.join(base_path, f) + gen = os.path.join(test_path, f) + if f != test["netcdf_output_file"]: with open(base, "r") as file1, open(gen, "r") as file2: # Skip first line next(file1) @@ -981,14 +978,14 @@ def test_disv01b(function_tmpdir, example_data_path): strt_longname = "starting head" ws = function_tmpdir / "ws" - for base_folder, test_info in tests.items(): - data_path = os.path.join(data_path_base, base_folder, test_info["base_sim_dir"]) + for dirname, test in tests.items(): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) # copy example data into working directory - base_model_folder = os.path.join(ws, f"{base_folder}_base") - test_model_folder = os.path.join(ws, f"{base_folder}_test") - shutil.copytree(data_path, base_model_folder) - os.mkdir(test_model_folder) + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) + os.mkdir(test_path) # create discretization disv = VertexGrid( @@ -1006,8 +1003,8 @@ def test_disv01b(function_tmpdir, example_data_path): ds = create_dataset( "gwf6", "disv01b", - test_info["netcdf_type"], - test_info["netcdf_output_file"], + test["netcdf_type"], + test["netcdf_output_file"], disv, ) @@ -1028,12 +1025,12 @@ def test_disv01b(function_tmpdir, example_data_path): ds.create_array("ic", "strt", strt, ["nlay", "ncpl"], strt_longname) # write to netcdf - ds.write(test_model_folder) + ds.write(test_path) # compare compare_netcdf( - os.path.join(base_model_folder, test_info["netcdf_output_file"]), - os.path.join(test_model_folder, test_info["netcdf_output_file"]), + os.path.join(base_path, test["netcdf_output_file"]), + os.path.join(test_path, test["netcdf_output_file"]), projection=True, ) @@ -1189,3 +1186,90 @@ def test_disv_transform(function_tmpdir, example_data_path): compare_netcdf_data(cmp_pth / f"tri.{nc_type}.nc", mf6_ws / "tri.in.nc") # compare_netcdf_data(cmp_pth / f"tri.{nc_type}.nc", vertex_ws / "tri.nc") + + +@pytest.mark.regression +def test_utlncf_load(function_tmpdir, example_data_path): + data_path_base = example_data_path / "mf6" / "netcdf" + tests = { + "test_utlncf_load": { + "base_sim_dir": "disv01b", + "netcdf_output_file": "disv01b.in.nc", + }, + } + ws = function_tmpdir / "ws" + for dirname, test in tests.items(): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) + + # copy example data into working directory + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) + + # load example + sim = flopy.mf6.MFSimulation.load(sim_ws=base_path) + + # set simulation path and write simulation + sim.set_sim_path(test_path) + sim.write_simulation(netcdf="mesh2d") + + # compare generated files + gen_files = [ + f + for f in os.listdir(test_path) + if os.path.isfile(os.path.join(test_path, f)) + ] + base_files = [ + f + for f in os.listdir(base_path) + if os.path.isfile(os.path.join(base_path, f)) + ] + + assert len(gen_files) == len(base_files) + for f in base_files: + base = os.path.join(base_path, f) + gen = os.path.join(test_path, f) + if f != test["netcdf_output_file"]: + with open(base, "r") as file1, open(gen, "r") as file2: + # Skip first line + next(file1) + next(file2) + + for line1, line2 in zip(file1, file2): + if line1.lower().startswith(" wkt"): + break + assert line1.lower() == line2.lower() + else: + compare_netcdf_data(base, gen) + + +@pytest.mark.regression +def test_utlncf_create(function_tmpdir, example_data_path): + data_path_base = example_data_path / "mf6" / "netcdf" + tests = { + "test_utlncf_create": { + "base_sim_dir": "disv01b", + "netcdf_output_file": "disv01b.in.nc", + }, + } + ws = function_tmpdir / "ws" + for dirname, test in tests.items(): + data_path = os.path.join(data_path_base, dirname, test["base_sim_dir"]) + + # copy example data into working directory + base_path = os.path.join(ws, f"{dirname}_base") + test_path = os.path.join(ws, f"{dirname}_test") + shutil.copytree(data_path, base_path) + + # load example + sim = flopy.mf6.MFSimulation.load(sim_ws=base_path) + + # set simulation path and write simulation + sim.set_sim_path(test_path) + sim.write_simulation(netcdf="mesh2d") + + # compare generated files + compare_path = os.path.join(base_path, "compare") + base = os.path.join(compare_path, "disv01b.in.nc") + gen = os.path.join(test_path, test["netcdf_output_file"]) + compare_netcdf_data(base, gen) diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/compare/disv01b.in.nc b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/compare/disv01b.in.nc new file mode 100644 index 000000000..ca87d5aed Binary files /dev/null and b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/compare/disv01b.in.nc differ diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.chd b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.chd new file mode 100644 index 000000000..bbf900dd2 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.chd @@ -0,0 +1,13 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options +END options + +BEGIN dimensions + MAXBOUND 2 +END dimensions + +BEGIN period 1 + 1 1 1.00000000E+00 + 1 9 0.00000000E+00 +END period 1 + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.disv b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.disv new file mode 100644 index 000000000..beab558e6 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.disv @@ -0,0 +1,63 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + EXPORT_ARRAY_NETCDF + NCF6 FILEIN disv01b.disv.ncf +END options + +BEGIN dimensions + NLAY 3 + NCPL 9 + NVERT 16 +END dimensions + +BEGIN griddata + top + INTERNAL FACTOR 1.0 + 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 + botm LAYERED + INTERNAL FACTOR 1.0 + -10.00000000 -10.00000000 -10.00000000 -10.00000000 -10.00000000 -10.00000000 -10.00000000 -10.00000000 -10.00000000 + INTERNAL FACTOR 1.0 + -20.00000000 -20.00000000 -20.00000000 -20.00000000 -20.00000000 -20.00000000 -20.00000000 -20.00000000 -20.00000000 + INTERNAL FACTOR 1.0 + -30.00000000 -30.00000000 -30.00000000 -30.00000000 -30.00000000 -30.00000000 -30.00000000 -30.00000000 -30.00000000 + idomain LAYERED + INTERNAL FACTOR 1 + 1 0 1 1 1 1 1 1 1 + INTERNAL FACTOR 1 + 1 1 1 1 1 1 1 1 1 + INTERNAL FACTOR 1 + 1 1 1 1 1 1 1 1 1 +END griddata + +BEGIN vertices + 1 1.00000000E+08 1.00000030E+08 + 2 1.00000010E+08 1.00000030E+08 + 3 1.00000020E+08 1.00000030E+08 + 4 1.00000030E+08 1.00000030E+08 + 5 1.00000000E+08 1.00000020E+08 + 6 1.00000010E+08 1.00000020E+08 + 7 1.00000020E+08 1.00000020E+08 + 8 1.00000030E+08 1.00000020E+08 + 9 1.00000000E+08 1.00000010E+08 + 10 1.00000010E+08 1.00000010E+08 + 11 1.00000020E+08 1.00000010E+08 + 12 1.00000030E+08 1.00000010E+08 + 13 1.00000000E+08 1.00000000E+08 + 14 1.00000010E+08 1.00000000E+08 + 15 1.00000020E+08 1.00000000E+08 + 16 1.00000030E+08 1.00000000E+08 +END vertices + +BEGIN cell2d + 1 1.00000005E+08 1.00000025E+08 4 1 2 6 5 + 2 1.00000015E+08 1.00000025E+08 4 2 3 7 6 + 3 1.00000025E+08 1.00000025E+08 4 3 4 8 7 + 4 1.00000005E+08 1.00000015E+08 4 5 6 10 9 + 5 1.00000015E+08 1.00000015E+08 4 6 7 11 10 + 6 1.00000025E+08 1.00000015E+08 4 7 8 12 11 + 7 1.00000005E+08 1.00000005E+08 4 9 10 14 13 + 8 1.00000015E+08 1.00000005E+08 4 10 11 15 14 + 9 1.00000025E+08 1.00000005E+08 4 11 12 16 15 +END cell2d + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.disv.ncf b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.disv.ncf new file mode 100644 index 000000000..49c0c6c8b --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.disv.ncf @@ -0,0 +1,9 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + wkt 'PROJCS["NAD83 / UTM zone 18N", GEOGCS["NAD83", DATUM["North_American_Datum_1983", SPHEROID["GRS 1980",6378137,298.257222101], TOWGS84[0,0,0,0,0,0,0]], PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4269"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-75], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1,AUTHORITY["EPSG","9001"]], AXIS["Easting",EAST], AXIS["Northing",NORTH], AUTHORITY["EPSG","26918"]]' + DEFLATE 9 + SHUFFLE + CHUNK_TIME 1 + CHUNK_FACE 3 +END options + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.ic b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.ic new file mode 100644 index 000000000..ca35b8f3e --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.ic @@ -0,0 +1,10 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + EXPORT_ARRAY_NETCDF +END options + +BEGIN griddata + strt + CONSTANT 0.00000000 +END griddata + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.ims b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.ims new file mode 100644 index 000000000..4b8778f35 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.ims @@ -0,0 +1,5 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + PRINT_OPTION summary +END options + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.nam b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.nam new file mode 100644 index 000000000..4ebb72535 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.nam @@ -0,0 +1,12 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options +END options + +BEGIN packages + DISV6 disv01b.disv disv + IC6 disv01b.ic ic + NPF6 disv01b.npf npf + CHD6 disv01b.chd chd_0 + OC6 disv01b.oc oc +END packages + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.npf b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.npf new file mode 100644 index 000000000..e0b4bc035 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.npf @@ -0,0 +1,12 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + EXPORT_ARRAY_NETCDF +END options + +BEGIN griddata + icelltype + CONSTANT 0 + k + CONSTANT 1.00000000 +END griddata + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.oc b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.oc new file mode 100644 index 000000000..10a8b2774 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.oc @@ -0,0 +1,9 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + HEAD FILEOUT disv01b.hds +END options + +BEGIN period 1 + SAVE HEAD ALL +END period 1 + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.tdis b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.tdis new file mode 100644 index 000000000..68234fa8b --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/disv01b.tdis @@ -0,0 +1,13 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + START_DATE_TIME 2041-01-01t00:00:00-05:00 +END options + +BEGIN dimensions + NPER 1 +END dimensions + +BEGIN perioddata + 1.00000000 1 1.00000000 +END perioddata + diff --git a/examples/data/mf6/netcdf/test_utlncf_create/disv01b/mfsim.nam b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/mfsim.nam new file mode 100644 index 000000000..b1e7500a5 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_create/disv01b/mfsim.nam @@ -0,0 +1,19 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options +END options + +BEGIN timing + TDIS6 disv01b.tdis +END timing + +BEGIN models + gwf6 disv01b.nam disv01b +END models + +BEGIN exchanges +END exchanges + +BEGIN solutiongroup 1 + ims6 disv01b.ims disv01b +END solutiongroup 1 + diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.chd b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.chd new file mode 100644 index 000000000..bbf900dd2 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.chd @@ -0,0 +1,13 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options +END options + +BEGIN dimensions + MAXBOUND 2 +END dimensions + +BEGIN period 1 + 1 1 1.00000000E+00 + 1 9 0.00000000E+00 +END period 1 + diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.disv b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.disv new file mode 100644 index 000000000..5cd9b5d63 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.disv @@ -0,0 +1,49 @@ +# NOT generated by Flopy +BEGIN options + EXPORT_ARRAY_NETCDF + NCF6 FILEIN disv01b.disv.ncf +END options + +BEGIN dimensions + NLAY 3 + NCPL 9 + NVERT 16 +END dimensions + +BEGIN griddata + top NETCDF + botm NETCDF + idomain NETCDF +END griddata + +BEGIN vertices + 1 1.00000000E+08 1.00000030E+08 + 2 1.00000010E+08 1.00000030E+08 + 3 1.00000020E+08 1.00000030E+08 + 4 1.00000030E+08 1.00000030E+08 + 5 1.00000000E+08 1.00000020E+08 + 6 1.00000010E+08 1.00000020E+08 + 7 1.00000020E+08 1.00000020E+08 + 8 1.00000030E+08 1.00000020E+08 + 9 1.00000000E+08 1.00000010E+08 + 10 1.00000010E+08 1.00000010E+08 + 11 1.00000020E+08 1.00000010E+08 + 12 1.00000030E+08 1.00000010E+08 + 13 1.00000000E+08 1.00000000E+08 + 14 1.00000010E+08 1.00000000E+08 + 15 1.00000020E+08 1.00000000E+08 + 16 1.00000030E+08 1.00000000E+08 +END vertices + +BEGIN cell2d + 1 1.00000005E+08 1.00000025E+08 4 1 2 6 5 + 2 1.00000015E+08 1.00000025E+08 4 2 3 7 6 + 3 1.00000025E+08 1.00000025E+08 4 3 4 8 7 + 4 1.00000005E+08 1.00000015E+08 4 5 6 10 9 + 5 1.00000015E+08 1.00000015E+08 4 6 7 11 10 + 6 1.00000025E+08 1.00000015E+08 4 7 8 12 11 + 7 1.00000005E+08 1.00000005E+08 4 9 10 14 13 + 8 1.00000015E+08 1.00000005E+08 4 10 11 15 14 + 9 1.00000025E+08 1.00000005E+08 4 11 12 16 15 +END cell2d + diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.disv.ncf b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.disv.ncf new file mode 100644 index 000000000..49c0c6c8b --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.disv.ncf @@ -0,0 +1,9 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + wkt 'PROJCS["NAD83 / UTM zone 18N", GEOGCS["NAD83", DATUM["North_American_Datum_1983", SPHEROID["GRS 1980",6378137,298.257222101], TOWGS84[0,0,0,0,0,0,0]], PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4269"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-75], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1,AUTHORITY["EPSG","9001"]], AXIS["Easting",EAST], AXIS["Northing",NORTH], AUTHORITY["EPSG","26918"]]' + DEFLATE 9 + SHUFFLE + CHUNK_TIME 1 + CHUNK_FACE 3 +END options + diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.ic b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.ic new file mode 100644 index 000000000..5ef059bd0 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.ic @@ -0,0 +1,8 @@ +# NOT generated by Flopy +BEGIN options + EXPORT_ARRAY_NETCDF +END options + +BEGIN griddata + strt NETCDF +END griddata diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.ims b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.ims new file mode 100644 index 000000000..4b8778f35 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.ims @@ -0,0 +1,5 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + PRINT_OPTION summary +END options + diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.in.nc b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.in.nc new file mode 100644 index 000000000..ca87d5aed Binary files /dev/null and b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.in.nc differ diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.nam b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.nam new file mode 100644 index 000000000..456a735e7 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.nam @@ -0,0 +1,12 @@ +# NOT generated by Flopy +BEGIN options + NETCDF FILEIN disv01b.in.nc +END options + +BEGIN packages + DISV6 disv01b.disv disv + IC6 disv01b.ic ic + NPF6 disv01b.npf npf + CHD6 disv01b.chd chd_0 + OC6 disv01b.oc oc +END packages diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.npf b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.npf new file mode 100644 index 000000000..381aa4640 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.npf @@ -0,0 +1,9 @@ +# NOT generated by Flopy +BEGIN options + EXPORT_ARRAY_NETCDF +END options + +BEGIN griddata + icelltype NETCDF + k NETCDF +END griddata diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.oc b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.oc new file mode 100644 index 000000000..10a8b2774 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.oc @@ -0,0 +1,9 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + HEAD FILEOUT disv01b.hds +END options + +BEGIN period 1 + SAVE HEAD ALL +END period 1 + diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.tdis b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.tdis new file mode 100644 index 000000000..68234fa8b --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/disv01b.tdis @@ -0,0 +1,13 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options + START_DATE_TIME 2041-01-01t00:00:00-05:00 +END options + +BEGIN dimensions + NPER 1 +END dimensions + +BEGIN perioddata + 1.00000000 1 1.00000000 +END perioddata + diff --git a/examples/data/mf6/netcdf/test_utlncf_load/disv01b/mfsim.nam b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/mfsim.nam new file mode 100644 index 000000000..b1e7500a5 --- /dev/null +++ b/examples/data/mf6/netcdf/test_utlncf_load/disv01b/mfsim.nam @@ -0,0 +1,19 @@ +# File generated by Flopy version 3.10.0.dev1 on 02/05/2025 at 13:05:06. +BEGIN options +END options + +BEGIN timing + TDIS6 disv01b.tdis +END timing + +BEGIN models + gwf6 disv01b.nam disv01b +END models + +BEGIN exchanges +END exchanges + +BEGIN solutiongroup 1 + ims6 disv01b.ims disv01b +END solutiongroup 1 + diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py index 385fdaecc..064c5b064 100644 --- a/flopy/mf6/mfmodel.py +++ b/flopy/mf6/mfmodel.py @@ -949,6 +949,8 @@ def load_base( break if dis_type and dis_type in dis_str: nc_fpth = os.path.join(instance.model_ws, nc_filerecord[0][0]) + # TODO: verify a partial dataset (not all griddata vars) don't + # break everything instance._nc_dataset = open_dataset(nc_fpth, dis_type=dis_str[dis_type]) else: message = ( @@ -1333,7 +1335,6 @@ def write( self, ext_file_action=ExtFileAction.copy_relative_paths, netcdf=None, - to_cdl=False, ): """ Writes out model's package files. @@ -1347,13 +1348,11 @@ def write( netcdf : str Create model NetCDF file, of type specified, in which to store package griddata. 'mesh2d' and 'structured' are supported types. - to_cdl : bool - Generate text version of netcdf file (debug feature) - """ # write netcdf file if netcdf or self._nc_dataset is not None: + kwargs = {} if self._nc_dataset is None: from ..utils.model_netcdf import create_dataset @@ -1370,10 +1369,21 @@ def write( # reset data storage and populate netcdf file for pp in self.packagelist: - pp._set_netcdf_storage(self._nc_dataset, create=True) + if pp.package_type == "ncf": + kwargs["shuffle"] = pp.shuffle.get_data() + kwargs["deflate"] = pp.deflate.get_data() + kwargs["chunk_time"] = pp.chunk_time.get_data() + kwargs["chunk_face"] = pp.chunk_face.get_data() + kwargs["chunk_x"] = pp.chunk_x.get_data() + kwargs["chunk_y"] = pp.chunk_y.get_data() + kwargs["chunk_z"] = pp.chunk_z.get_data() + kwargs["wkt"] = pp.wkt.get_data() + pp._set_netcdf_storage( + self._nc_dataset, create=True + ) # write the dataset to netcdf - self._nc_dataset.write(self.model_ws) + self._nc_dataset.write(self.model_ws, **kwargs) # write name file if ( diff --git a/flopy/mf6/mfsimbase.py b/flopy/mf6/mfsimbase.py index 4f70c579e..475d12c87 100644 --- a/flopy/mf6/mfsimbase.py +++ b/flopy/mf6/mfsimbase.py @@ -1668,7 +1668,6 @@ def write_simulation( ext_file_action=ExtFileAction.copy_relative_paths, silent=False, netcdf=None, - to_cdl=False, ): """ Write the simulation to files. @@ -1684,9 +1683,6 @@ def write_simulation( netcdf : str Create model NetCDF file, of type specified, in which to store package griddata. 'mesh2d' and 'structured' are supported types. - to_cdl : bool - Generate text version of netcdf file (debug feature) - """ sim_data = self.simulation_data if not sim_data.max_columns_user_set: @@ -1756,7 +1752,6 @@ def write_simulation( model.write( ext_file_action=ext_file_action, netcdf=netcdf, - to_cdl=to_cdl, ) self.simulation_data.mfpath.set_last_accessed_path() diff --git a/flopy/utils/model_netcdf.py b/flopy/utils/model_netcdf.py index f975ea894..cda4a3a99 100644 --- a/flopy/utils/model_netcdf.py +++ b/flopy/utils/model_netcdf.py @@ -28,8 +28,7 @@ class ModelNetCDFDataset: Newly created datasets define coordinate or mesh variables corresponding to the type of NetCDF file specified. When the discretization crs attribute is valid, a projection - variable is also added to the dataset and candidate data - variables are associated with the grid. + variable is also added to the dataset. Data will be associated with the grid and projection when the relevant interfaces, e.g. create_array() and write(), @@ -38,9 +37,10 @@ class ModelNetCDFDataset: Additionally, these files can can be used as MODFLOW 6 model inputs for variables that define internal attributes - designated for that purpose, specifically "modflow6_input". - These attributes are managed internally when the supported - data interfaces are used. + designated for that purpose, specifically "modflow6_input" + and "modflow6_layer". These attributes are managed internally + for MODFLOW 6 models when the supported data interfaces are + used. """ def __init__(self): @@ -129,7 +129,7 @@ def create_array( Create a new array. Override this function in a derived class. Args: - package (str): A name of a data grouping in the file, typically + package (str): The name of a data grouping in the file, typically a package. Must be distinct for each grouping. If this dataset is associated with a modflow 6 model and this is a base package (dis, disv, npf, ic, etc.), use that name. If this is a stress @@ -181,17 +181,25 @@ def array(self, package: str, param: str, layer=None): """ raise NotImplementedError("array not implemented in base class") - def write(self, path: str) -> None: + def write(self, path: str, **kwargs) -> None: """ Write the data set to a NetCDF file. Args: path (str): A directory in which to write the file. + kwargs (dict): A dictionay of supported encodings to + apply to managed grid associated arrays. """ self._set_projection() self._set_modflow_attrs() nc_fpath = Path(path) / self._fname - self._dataset.to_netcdf(nc_fpath, format="NETCDF4", engine="netcdf4") + encoding = self._set_encoding(**kwargs) + if encoding is not None: + self._dataset.to_netcdf( + nc_fpath, format="NETCDF4", engine="netcdf4", encoding=encoding + ) + else: + self._dataset.to_netcdf(nc_fpath, format="NETCDF4", engine="netcdf4") def path(self, package: str, param: str): return f"{self._modelname.lower()}/{package.lower()}/{param.lower()}" @@ -206,7 +214,7 @@ def _set_grid(self, dis): """ raise NotImplementedError("_set_grid not implemented in base class") - def _set_coords(self, dis): + def _set_coords(self, crs): """ Define coordinate variables associated with the NetCDF file type. Override this function in a derived class. @@ -216,23 +224,110 @@ def _set_coords(self, dis): """ raise NotImplementedError("_set_coords not implemented in base class") + def _set_encoding(self, **kwargs): + if self._dis is None: + return None + # encodings: { + # 'szip_coding', + # 'shuffle', + # 'fletcher32', + # 'quantize_mode', + # 'least_significant_digit', + # 'endian', + # 'szip_pixels_per_block', + # '_FillValue', + # 'compression', + # 'blosc_shuffle', + # 'zlib', + # 'significant_digits', + # 'complevel', + # 'dtype', + # 'contiguous', + # 'chunksizes'} + encoding = {} + encodes = {} + + deflate = kwargs.pop("deflate", None) + shuffle = kwargs.pop("shuffle", None) + chunk_time = kwargs.pop("chunk_time", None) + chunk_face = kwargs.pop("chunk_face", None) + chunk_x = kwargs.pop("chunk_x", None) + chunk_y = kwargs.pop("chunk_y", None) + chunk_z = kwargs.pop("chunk_z", None) + + if deflate: + encodes["zlib"] = True + encodes["complevel"] = deflate + if shuffle: + encodes["shuffle"] = True + + for path in self._tags: + for l in self._tags[path]: + if chunk_face and self._nc_type == "mesh2d": + codes = dict(encodes) + dims = self._dataset[self._tags[path][l]].dims + if "nmesh_face" in dims: + codes["chunksizes"] = [chunk_face] + encoding[self._tags[path][l]] = codes + elif self._nc_type == "structured" and chunk_x and chunk_y and chunk_z: + codes = dict(encodes) + dims = self._dataset[self._tags[path][l]].dims + if "x" in dims and "y" in dims: + if "z" in dims: + codes["chunksizes"] = [chunk_z, chunk_y, chunk_x] + else: + codes["chunksizes"] = [chunk_y, chunk_x] + encoding[self._tags[path][l]] = codes + + else: + encoding[self._tags[path][l]] = encodes + + return encoding + def _set_projection(self): - if self._dis and self._dis.crs: - crs = CRS.from_user_input(self._dis.crs) + if self._dis: + crs = None + wkt = None + if self._dis.crs: + try: + crs = CRS.from_user_input(self._dis.crs) + wkt = crs.to_wkt() + except Exception as e: + # crs = None + # wkt = None + warnings.warn( + f"Cannot generate CRS object from user input: {e}", + UserWarning, + ) + + if wkt is not None: + # add projection variable + self._dataset = self._dataset.assign( + {"projection": ([], np.int32(1))} + ) + if self._nc_type == "structured": + self._dataset["projection"].attrs["crs_wkt"] = wkt + self._dataset["x"].attrs["grid_mapping"] = "projection" + self._dataset["y"].attrs["grid_mapping"] = "projection" + elif self._nc_type == "mesh2d": + self._dataset["projection"].attrs["wkt"] = wkt + self._dataset["mesh_node_x"].attrs["grid_mapping"] = ( + "projection" + ) + self._dataset["mesh_node_y"].attrs["grid_mapping"] = ( + "projection" + ) + self._dataset["mesh_face_x"].attrs["grid_mapping"] = ( + "projection" + ) + self._dataset["mesh_face_y"].attrs["grid_mapping"] = ( + "projection" + ) + + # update coords based on crs coords = self._set_coords(crs) - wkt = crs.to_wkt() - self._dataset = self._dataset.assign({"projection": ([], np.int32(1))}) - if self._nc_type == "structured": - self._dataset["projection"].attrs["crs_wkt"] = wkt - self._dataset["x"].attrs["grid_mapping"] = "projection" - self._dataset["y"].attrs["grid_mapping"] = "projection" - elif self._nc_type == "mesh2d": - self._dataset["projection"].attrs["wkt"] = wkt - self._dataset["mesh_node_x"].attrs["grid_mapping"] = "projection" - self._dataset["mesh_node_y"].attrs["grid_mapping"] = "projection" - self._dataset["mesh_face_x"].attrs["grid_mapping"] = "projection" - self._dataset["mesh_face_y"].attrs["grid_mapping"] = "projection" + # set grid_mapping and coordinates attributes for p in self._tags: for l in self._tags[p]: dims = self._dataset[self._tags[p][l]].dims @@ -240,9 +335,10 @@ def _set_projection(self): self._nc_type == "mesh2d" and ("nmesh_face" in dims or "nmesh_node" in dims) ): - self._dataset[self._tags[p][l]].attrs["grid_mapping"] = ( - "projection" - ) + if crs is not None: + self._dataset[self._tags[p][l]].attrs["grid_mapping"] = ( + "projection" + ) if coords is not None: self._dataset[self._tags[p][l]].attrs["coordinates"] = ( coords @@ -279,10 +375,6 @@ def _set_mapping(self): if "modflow6_input" in da.attrs: path = da.attrs["modflow6_input"].lower() - tokens = path.split("/") - assert len(tokens) == 3 - assert tokens[0] == self._modelname - if "modflow6_layer" in da.attrs: layer = da.attrs["modflow6_layer"] # convert indexing to 0-based @@ -388,7 +480,6 @@ def _create_layered_array( else: ln = longname self._dataset[layer_vname].attrs["long_name"] = ln - self._dataset[layer_vname].attrs["coordinates"] = "mesh_face_x mesh_face_y" if path not in self._tags: self._tags[path] = {} if layer in self._tags[path]: @@ -479,7 +570,7 @@ def _set_grid(self, dis): # y = dis.ycellcenters[:, 0] x = dis.xoffset + dis.xycenters[0] y = dis.yoffset + dis.xycenters[1] - z = list(range(1, dis.nlay + 1)) + z = [float(x) for x in range(1, dis.nlay + 1)] # create coordinate vars var_d = {"z": (["z"], z), "y": (["y"], y), "x": (["x"], x)} @@ -504,6 +595,9 @@ def _set_grid(self, dis): self._dataset["x"].attrs["bounds"] = "x_bnds" def _set_coords(self, crs): + if crs is None: + return "x y" + lats = [] lons = [] xdim = self._dataset.dims["x"] @@ -714,8 +808,8 @@ def _set_grid(self, dis): self._dataset["mesh_face_nodes"].attrs["_FillValue"] = FILLNA_INT32 self._dataset["mesh_face_nodes"].attrs["start_index"] = np.int32(1) - def _set_coords(self, proj): - return None + def _set_coords(self, crs): + return "mesh_face_x mesh_face_y" class DisvNetCDFMesh2d(ModelNetCDFDataset): @@ -857,8 +951,8 @@ def _set_grid(self, disv): self._dataset["mesh_face_nodes"].attrs["_FillValue"] = FILLNA_INT32 self._dataset["mesh_face_nodes"].attrs["start_index"] = np.int32(1) - def _set_coords(self, proj): - return None + def _set_coords(self, crs): + return "mesh_face_x mesh_face_y" def open_dataset(nc_fpth: str, dis_type: str) -> ModelNetCDFDataset: @@ -905,7 +999,6 @@ def open_dataset(nc_fpth: str, dis_type: str) -> ModelNetCDFDataset: dataset.close() if nc_dataset: - fpth = Path(nc_fpth).resolve() nc_dataset.open(fpth) else: raise Exception(