docs: freyberg mfusg & gsf file example

* add utility function to read unstructured grid props from gsf file
* add example notebook demonstrating freyberg mfusg model from gsf

autotest/test_grid.py

@@ -12,6 +12,7 @@
 from pytest_cases import parametrize_with_cases
 
 from flopy.discretization import StructuredGrid, UnstructuredGrid, VertexGrid
+from flopy.discretization.unstructuredgrid import unstructured_props_from_gsf
 from flopy.mf6 import MFSimulation
 from flopy.modflow import Modflow, ModflowDis
 from flopy.utils.cvfdutil import gridlist_to_disv_gridprops, to_cvfd
@@ -910,3 +911,25 @@ def test_get_lni_unstructured(grid):
             )
         )
         assert csum[layer] + i == nn
+
+
+def test_unstructured_props_from_gsf(example_data_path):
+    model_path = example_data_path / "freyberg_usg"
+    verts, iverts, xc, yc, ncpl, top, botm = unstructured_props_from_gsf(
+        str(model_path / "freyberg.usg.gsf")
+    )
+
+    assert len(verts) == 7020
+    assert verts[0] == [0, 619741.4, 3373170.9]
+    assert verts[-1] == [7019, 624741.4, 3363170.9]
+    assert len(iverts) == 4497
+    assert iverts[0] == [0, 36, 37, 1, 0]
+    assert iverts[-1] == [5227, 5263, 5264, 5228, 5227]
+    assert np.max(xc) == 624616.4
+    assert np.max(yc) == 3373139.65
+    assert np.array_equal(ncpl, [1499, 1499, 1499])
+    assert xc.shape == yc.shape == top.shape == botm.shape == (4497,)
+    assert np.min(top) == -1049.99
+    assert np.max(top) == 206871.7
+    assert np.min(botm) == -508.745
+    assert np.max(botm) == 33.31749