tf_hidrogeo

{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":6702999,"sourceType":"datasetVersion","datasetId":3863347}],"dockerImageVersionId":30626,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"source":"<a href=\"https://www.kaggle.com/code/jorgeluiscruzquispe/tf-hidrogeologia-t?scriptVersionId=171431909\" target=\"_blank\"><img align=\"left\" alt=\"Kaggle\" title=\"Open in Kaggle\" src=\"https://kaggle.com/static/images/open-in-kaggle.svg\"></a>","metadata":{},"cell_type":"markdown"},{"cell_type":"markdown","source":"# <center>Ejemplo de geoprocesamiento de cuencas<center>\n\n---\n\nFrom Hughes, Joseph D., Langevin, Christian D., Paulinski, Scott R., Larsen, Joshua D., and Brakenhoff, David, 2023, Flujos de trabajo FloPy para crear modelos MODFLOW estructurados y no estructurados: aguas subterráneas, https://doi.org/10.1111/gwat.13327","metadata":{}},{"cell_type":"raw","source":"# !pip install git+https://github.com/modflowpy/flopy.git\n!pip install -q flopy\n!get-modflow :flopy","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:18.162925Z","iopub.execute_input":"2024-04-11T01:15:18.163971Z","iopub.status.idle":"2024-04-11T01:15:38.003947Z","shell.execute_reply.started":"2024-04-11T01:15:18.163923Z","shell.execute_reply":"2024-04-11T01:15:38.002585Z"}}},{"cell_type":"code","source":"import os\nimport pathlib as pl\nimport sys","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.status.busy":"2024-04-11T01:15:38.006593Z","iopub.execute_input":"2024-04-11T01:15:38.007157Z","iopub.status.idle":"2024-04-11T01:15:38.014618Z","shell.execute_reply.started":"2024-04-11T01:15:38.007103Z","shell.execute_reply":"2024-04-11T01:15:38.013141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib as mpl\nimport matplotlib.gridspec as gridspec\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom shapely.geometry import LineString","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.016992Z","iopub.execute_input":"2024-04-11T01:15:38.017425Z","iopub.status.idle":"2024-04-11T01:15:38.568141Z","shell.execute_reply.started":"2024-04-11T01:15:38.017388Z","shell.execute_reply":"2024-04-11T01:15:38.567069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import flopy\nimport flopy.plot.styles as styles\nfrom flopy.discretization import StructuredGrid, VertexGrid\nfrom flopy.utils.gridgen import Gridgen\nfrom flopy.utils.triangle import Triangle\nfrom flopy.utils.voronoi import VoronoiGrid","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.569912Z","iopub.execute_input":"2024-04-11T01:15:38.571093Z","iopub.status.idle":"2024-04-11T01:15:38.585373Z","shell.execute_reply.started":"2024-04-11T01:15:38.571038Z","shell.execute_reply":"2024-04-11T01:15:38.58431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(sys.version)\nprint(f\"numpy version: {np.__version__}\")\nprint(f\"matplotlib version: {mpl.__version__}\")\nprint(f\"flopy version: {flopy.__version__}\")","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.5869Z","iopub.execute_input":"2024-04-11T01:15:38.588053Z","iopub.status.idle":"2024-04-11T01:15:38.598106Z","shell.execute_reply.started":"2024-04-11T01:15:38.588012Z","shell.execute_reply":"2024-04-11T01:15:38.597073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# importar toda la información de estilo de trazado desde defaults.py\nsys.path.append(\"../common\")\nimport numpy as np\nimport shapely\nfrom shapely.geometry import Polygon\n\nfrom flopy.utils.gridintersect import GridIntersect\n\ngeometries = {\n    \"boundary\": \"\"\"1.868012422360248456e+05 4.695652173913043953e+04\n    1.790372670807453396e+05 5.204968944099379587e+04\n    1.729813664596273447e+05 5.590062111801243009e+04\n    1.672360248447204940e+05 5.987577639751553215e+04\n    1.631987577639751253e+05 6.335403726708075556e+04\n    1.563664596273291972e+05 6.819875776397516893e+04\n    1.509316770186335489e+05 7.229813664596274612e+04\n    1.453416149068323139e+05 7.527950310559007630e+04\n    1.395962732919254631e+05 7.627329192546584818e+04\n    1.357142857142857101e+05 7.664596273291927355e+04\n    1.329192546583850926e+05 7.751552795031057030e+04\n    1.268633540372670832e+05 8.062111801242237561e+04\n    1.218944099378881947e+05 8.285714285714286962e+04\n    1.145962732919254486e+05 8.571428571428572468e+04\n    1.069875776397515583e+05 8.869565217391305487e+04\n    1.023291925465838431e+05 8.931677018633540138e+04\n    9.456521739130433707e+04 9.068322981366459862e+04\n    8.804347826086955320e+04 9.080745341614908830e+04\n    7.950310559006211406e+04 9.267080745341615693e+04\n    7.562111801242236106e+04 9.391304347826087906e+04\n    6.692546583850930620e+04 9.602484472049689793e+04\n    5.667701863354037778e+04 9.763975155279504543e+04\n    4.906832298136646568e+04 9.689440993788820924e+04\n    3.897515527950309479e+04 9.540372670807455142e+04\n    3.167701863354036323e+04 9.304347826086958230e+04\n    2.375776397515527788e+04 8.757763975155279331e+04\n    1.847826086956521613e+04 8.161490683229814749e+04\n    1.164596273291925172e+04 7.739130434782608063e+04\n    6.211180124223596977e+03 7.055900621118013805e+04\n    4.347826086956512881e+03 6.422360248447205959e+04\n    1.863354037267072272e+03 6.037267080745341809e+04\n    2.639751552795024509e+03 5.602484472049689793e+04\n    1.552795031055893560e+03 5.279503105590062478e+04\n    7.763975155279410956e+02 4.186335403726709046e+04\n    2.018633540372667312e+03 3.813664596273292409e+04\n    6.055900621118013078e+03 3.341614906832297856e+04\n    1.335403726708074100e+04 2.782608695652173992e+04\n    2.577639751552794405e+04 2.086956521739130767e+04\n    3.416149068322980747e+04 1.763975155279503815e+04\n    4.642857142857142753e+04 1.440993788819875044e+04\n    5.636645962732918997e+04 1.130434782608694877e+04\n    6.459627329192546313e+04 9.813664596273290954e+03\n    8.555900621118012350e+04 6.832298136645956220e+03\n    9.829192546583850344e+04 5.093167701863346338e+03\n    1.085403726708074391e+05 4.347826086956525614e+03\n    1.200310559006211115e+05 4.223602484472040487e+03\n    1.296583850931677007e+05 4.347826086956525614e+03\n    1.354037267080745369e+05 5.590062111801232277e+03\n    1.467391304347825935e+05 1.267080745341615875e+04\n    1.563664596273291972e+05 1.937888198757762802e+04\n    1.630434782608695677e+05 2.198757763975155467e+04\n    1.694099378881987650e+05 2.434782608695652743e+04\n    1.782608695652173774e+05 2.981366459627329095e+04\n    1.833850931677018234e+05 3.180124223602484562e+04\n    1.868012422360248456e+05 3.577639751552795497e+04\"\"\",\n    \"streamseg1\": \"\"\"1.868012422360248456e+05 4.086956521739130403e+04\n    1.824534161490683327e+05 4.086956521739130403e+04\n    1.770186335403726553e+05 4.124223602484472940e+04\n    1.737577639751552779e+05 4.186335403726709046e+04\n    1.703416149068323139e+05 4.310559006211180531e+04\n    1.670807453416148783e+05 4.397515527950310934e+04\n    1.636645962732919143e+05 4.484472049689441337e+04\n    1.590062111801242281e+05 4.559006211180124228e+04\n    1.555900621118012350e+05 4.559006211180124228e+04\n    1.510869565217391064e+05 4.546583850931677443e+04\n    1.479813664596273156e+05 4.534161490683229931e+04\n    1.453416149068323139e+05 4.496894409937888850e+04\n    1.377329192546583654e+05 4.447204968944099528e+04\n    1.326086956521739194e+05 4.447204968944099528e+04\n    1.285714285714285652e+05 4.434782608695652743e+04\n    1.245341614906832110e+05 4.472049689440993825e+04\n    1.215838509316770069e+05 4.509316770186335634e+04\n    1.161490683229813585e+05 4.509316770186335634e+04\n    1.125776397515527933e+05 4.459627329192547040e+04\n    1.074534161490683036e+05 4.385093167701864149e+04\n    1.018633540372670686e+05 4.347826086956522340e+04\n    9.798136645962731563e+04 4.360248447204969125e+04\n    9.223602484472049400e+04 4.310559006211180531e+04\n    8.602484472049689793e+04 4.198757763975155831e+04\n    7.981366459627327276e+04 4.173913043478261534e+04\n    7.468944099378881219e+04 4.248447204968944425e+04\n    7.034161490683228476e+04 4.385093167701864149e+04\n    6.785714285714285506e+04 4.621118012422360334e+04\n    6.583850931677018525e+04 4.919254658385094081e+04\n    6.319875776397513982e+04 5.192546583850932075e+04\n    6.009316770186335634e+04 5.677018633540373412e+04\n    5.605590062111800216e+04 5.950310559006211406e+04\n    5.279503105590060295e+04 6.124223602484472940e+04\n    4.751552795031056303e+04 6.211180124223603343e+04\n    3.990683229813664366e+04 6.335403726708075556e+04\n    3.276397515527949508e+04 6.409937888198757719e+04\n    2.934782608695651652e+04 6.509316770186336362e+04\n    2.546583850931676716e+04 6.832298136645962950e+04\"\"\",\n    \"streamseg2\": \"\"\"7.025161490683228476e+04 4.375093167701864149e+04\n    6.816770186335404287e+04 4.273291925465839449e+04\n    6.490683229813665093e+04 4.211180124223603343e+04\n    6.164596273291925900e+04 4.173913043478262261e+04\n    5.776397515527951327e+04 4.124223602484472940e+04\n    5.450310559006211406e+04 4.049689440993789322e+04\n    4.984472049689442065e+04 3.937888198757764621e+04\n    4.534161490683231386e+04 3.801242236024845624e+04\n    4.114906832298137306e+04 3.664596273291926627e+04\n    3.913043478260868869e+04 3.565217391304348712e+04\n    3.649068322981366509e+04 3.416149068322981475e+04\n    3.322981366459628043e+04 3.242236024844721760e+04\n    3.012422360248447148e+04 3.105590062111801672e+04\n    2.608695652173913550e+04 2.957521739130435890e+04\"\"\",\n    \"streamseg3\": \"\"\"1.059006211180124228e+05 4.335403726708074828e+04\n    1.029503105590062187e+05 4.223602484472050128e+04\n    1.004658385093167890e+05 4.024844720496894297e+04\n    9.937888198757765349e+04 3.788819875776398112e+04\n    9.627329192546584818e+04 3.490683229813664366e+04\n    9.285714285714286962e+04 3.316770186335403559e+04\n    8.897515527950311662e+04 3.093167701863354159e+04\n    8.338509316770188161e+04 2.795031055900621504e+04\n    7.872670807453416637e+04 2.670807453416148928e+04\n    7.329192546583851799e+04 2.385093167701863058e+04\n    6.863354037267081731e+04 2.111801242236025064e+04\n    6.304347826086958230e+04 1.863354037267081003e+04\"\"\",\n    \"streamseg4\": \"\"\"1.371118012422360480e+05 4.472049689440994553e+04\n    1.321428571428571595e+05 4.720496894409938250e+04\n    1.285714285714285652e+05 4.981366459627330187e+04\n    1.243788819875776535e+05 5.341614906832298584e+04\n    1.189440993788819906e+05 5.540372670807454415e+04\n    1.125776397515527933e+05 5.627329192546584818e+04\n    1.065217391304347839e+05 5.726708074534162733e+04\n    1.020186335403726698e+05 5.913043478260870324e+04\n    9.409937888198759174e+04 6.273291925465840177e+04\n    9.192546583850932075e+04 6.633540372670808574e+04\n    8.881987577639751544e+04 7.242236024844722124e+04\n    8.586956521739131131e+04 7.552795031055902655e+04\n    8.369565217391305487e+04 7.962732919254660374e+04\"\"\",\n}\n\ndef string2geom(geostring, conversion=None):\n    if conversion is None:\n        multiplier = 1.0\n    else:\n        multiplier = float(conversion)\n    res = []\n    for line in geostring.split(\"\\n\"):\n        line = line.strip()\n        line = line.split(\" \")\n        x = float(line[0]) * multiplier\n        y = float(line[1]) * multiplier\n        res.append((x, y))\n    return res\n\ndef densify_geometry(line, step, keep_internal_nodes=True):\n    xy = []  # list of tuple of coordinates\n    lines_strings = []\n    if keep_internal_nodes:\n        for idx in range(1, len(line)):\n            lines_strings.append(\n                shapely.geometry.LineString(line[idx - 1 : idx + 1])\n            )\n    else:\n        lines_strings = [shapely.geometry.LineString(line)]\n\n    for line_string in lines_strings:\n        length_m = line_string.length  # get the length\n        for distance in np.arange(0, length_m + step, step):\n            point = line_string.interpolate(distance)\n            xy_tuple = (point.x, point.y)\n            if xy_tuple not in xy:\n                xy.append(xy_tuple)\n        # make sure the end point is in xy\n        if keep_internal_nodes:\n            xy_tuple = line_string.coords[-1]\n            if xy_tuple not in xy:\n                xy.append(xy_tuple)\n    return xy\n\n# Función para configurar el área del modelo activo e inactivo.\ndef set_idomain(grid, boundary):\n    ix = GridIntersect(grid, method=\"vertex\", rtree=True)\n    result = ix.intersect(Polygon(boundary))\n    idx = [coords for coords in result.cellids]\n    idx = np.array(idx, dtype=int)\n    nr = idx.shape[0]\n    if idx.ndim == 1:\n        idx = idx.reshape((nr, 1))\n    idx = tuple([idx[:, i] for i in range(idx.shape[1])])\n    idomain = np.zeros(grid.shape[1:], dtype=int)\n    idomain[idx] = 1\n    idomain = idomain.reshape(grid.shape)\n    grid.idomain = idomain","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.600074Z","iopub.execute_input":"2024-04-11T01:15:38.600844Z","iopub.status.idle":"2024-04-11T01:15:38.629384Z","shell.execute_reply.started":"2024-04-11T01:15:38.600803Z","shell.execute_reply":"2024-04-11T01:15:38.627767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"tamaño de figura basica","metadata":{}},{"cell_type":"code","source":"figwidth = 360  # 90 # mm\nfigwidth = figwidth / 10 / 2.54  # inches\nfigheight = figwidth\nfigsize = (figwidth, figheight)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.631535Z","iopub.execute_input":"2024-04-11T01:15:38.632549Z","iopub.status.idle":"2024-04-11T01:15:38.648086Z","shell.execute_reply.started":"2024-04-11T01:15:38.632505Z","shell.execute_reply":"2024-04-11T01:15:38.646977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"configuración de las figuras","metadata":{}},{"cell_type":"code","source":"levels = np.arange(10, 110, 10)\ncontour_color = \"black\"\ncontour_style = \"--\"\ncontour_dict = {\n    \"levels\": levels,\n    \"linewidths\": 0.5,\n    \"colors\": contour_color,\n    \"linestyles\": contour_style,\n}\nclabel_dict = {\n    \"inline\": True,\n    \"fmt\": \"%1.0f\",\n    \"fontsize\": 6,\n    \"inline_spacing\": 0.5,\n}\nfont_dict = {\n    \"fontsize\": 5,\n    \"color\": \"black\",\n}\ngrid_dict = {\n    \"lw\": 0.25,\n    \"color\": \"0.5\",\n}\nrefinement_dict = {\n    \"color\": \"magenta\",\n    \"ls\": \":\",\n    \"lw\": 1.0,\n}\nriver_dict = {\n    \"color\": \"blue\",\n    \"linestyle\": \"-\",\n    \"linewidth\": 1,\n}\nintersection_cmap = \"Pastel2\"\ntemp_cmap = mpl.colormaps[intersection_cmap]\nintersection_rgba = temp_cmap(0)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.649917Z","iopub.execute_input":"2024-04-11T01:15:38.650645Z","iopub.status.idle":"2024-04-11T01:15:38.667466Z","shell.execute_reply.started":"2024-04-11T01:15:38.650602Z","shell.execute_reply":"2024-04-11T01:15:38.666212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Se crea una carpeta temporal para los outputs","metadata":{}},{"cell_type":"code","source":"temp_path = \"./temp\"\nif not os.path.isdir(temp_path):\n    os.mkdir(temp_path)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.66976Z","iopub.execute_input":"2024-04-11T01:15:38.670349Z","iopub.status.idle":"2024-04-11T01:15:38.684039Z","shell.execute_reply.started":"2024-04-11T01:15:38.670293Z","shell.execute_reply":"2024-04-11T01:15:38.682484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cargando raster de la cuenca","metadata":{}},{"cell_type":"code","source":"dem_ascii = \"/kaggle/input/ascii-file/fine_topo.asc\"\nfine_topo = flopy.utils.Raster.load(dem_ascii)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.686237Z","iopub.execute_input":"2024-04-11T01:15:38.687167Z","iopub.status.idle":"2024-04-11T01:15:38.717489Z","shell.execute_reply.started":"2024-04-11T01:15:38.68709Z","shell.execute_reply":"2024-04-11T01:15:38.716093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Definir el tamaño y la extensión del problema.","metadata":{}},{"cell_type":"code","source":"Lx = 180000\nLy = 100000\nextent = (0, Lx, 0, Ly)\nvmin, vmax = 0.0, 100.0","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.719034Z","iopub.execute_input":"2024-04-11T01:15:38.719393Z","iopub.status.idle":"2024-04-11T01:15:38.725228Z","shell.execute_reply.started":"2024-04-11T01:15:38.719363Z","shell.execute_reply":"2024-04-11T01:15:38.724074Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cree datos de límites y ríos que se utilizarán para refinar las cuadrículas no estructuradas.","metadata":{}},{"cell_type":"code","source":"boundary_polygon = string2geom(geometries[\"boundary\"])\nbp = np.array(boundary_polygon)\n\nsgs = [string2geom(geometries[f\"streamseg{i}\"]) for i in range(1, 5)]\n\nfig = plt.figure(figsize=figsize)\nax = fig.add_subplot()\nax.set_aspect(\"equal\")\n\nriv_colors = (\"blue\", \"cyan\", \"green\", \"orange\", \"red\")\n\nax.plot(bp[:, 0], bp[:, 1], \"ko-\")\nfor idx, sg in enumerate(sgs):\n    sa = np.array(sg)\n    ax.plot(sa[:, 0], sa[:, 1], color=riv_colors[idx], lw=0.75, marker=\"o\")","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:38.726712Z","iopub.execute_input":"2024-04-11T01:15:38.727108Z","iopub.status.idle":"2024-04-11T01:15:39.103434Z","shell.execute_reply.started":"2024-04-11T01:15:38.727075Z","shell.execute_reply":"2024-04-11T01:15:39.102308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Cuadrículas estructuradas\n#### Cuadrícula estructurada con espaciado constante entre filas y columnas.","metadata":{}},{"cell_type":"code","source":"dx = dy = 1666.66666667\n# dx = dy = 500\nnlay = 1\nnrow = int(Ly / dy) + 1\nncol = int(Lx / dx) + 1\ndelr = np.array(ncol * [dx])\ndelc = np.array(nrow * [dy])\ntop = np.ones((nrow, ncol)) * 1000.0\nbotm = np.ones((nlay, nrow, ncol)) * -100.0\n\nstruct_grid = StructuredGrid(nlay=nlay, delr=delr, delc=delc, xoff=0.0, yoff=0.0, top=top, botm=botm)\nset_idomain(struct_grid, boundary_polygon)","metadata":{"tags":[],"execution":{"iopub.status.busy":"2024-04-11T01:15:39.105097Z","iopub.execute_input":"2024-04-11T01:15:39.105453Z","iopub.status.idle":"2024-04-11T01:15:40.626404Z","shell.execute_reply.started":"2024-04-11T01:15:39.105422Z","shell.execute_reply":"2024-04-11T01:15:40.625141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_sg = fine_topo.resample_to_grid(\n    struct_grid,\n    band=fine_topo.bands[0],\n    method=\"linear\",\n    extrapolate_edges=True,)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:40.628241Z","iopub.execute_input":"2024-04-11T01:15:40.628697Z","iopub.status.idle":"2024-04-11T01:15:41.217209Z","shell.execute_reply.started":"2024-04-11T01:15:40.628653Z","shell.execute_reply":"2024-04-11T01:15:41.215677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ixs = flopy.utils.GridIntersect(struct_grid, method=\"structured\")\ncellids = []\nfor sg in sgs:\n    v = ixs.intersect(LineString(sg), sort_by_cellid=True)\n    cellids += v[\"cellids\"].tolist()\nintersection_sg = np.zeros(struct_grid.shape[1:])\nfor loc in cellids:\n    intersection_sg[loc] = 1","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:41.219088Z","iopub.execute_input":"2024-04-11T01:15:41.224835Z","iopub.status.idle":"2024-04-11T01:15:41.387485Z","shell.execute_reply.started":"2024-04-11T01:15:41.224761Z","shell.execute_reply":"2024-04-11T01:15:41.386174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=figsize)\nax = fig.add_subplot()\npmv = flopy.plot.PlotMapView(modelgrid=struct_grid)\nax.set_aspect(\"equal\")\npmv.plot_array(top_sg)\npmv.plot_array(\n    intersection_sg,\n    masked_values=[\n        0,\n    ],\n    alpha=0.2,\n    cmap=\"Reds_r\",\n)\npmv.plot_grid(lw=0.25, color=\"0.5\")\ncg = pmv.contour_array(top_sg, levels=levels, linewidths=0.3, colors=\"0.75\")\npmv.plot_inactive()\n\nax.plot(bp[:, 0], bp[:, 1], \"k-\")\nfor sg in sgs:\n    sa = np.array(sg)\n    ax.plot(sa[:, 0], sa[:, 1], \"b-\")","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:41.389207Z","iopub.execute_input":"2024-04-11T01:15:41.389576Z","iopub.status.idle":"2024-04-11T01:15:41.924927Z","shell.execute_reply.started":"2024-04-11T01:15:41.389544Z","shell.execute_reply":"2024-04-11T01:15:41.923755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize = (15, 5))\nxs = flopy.plot.PlotCrossSection(modelgrid=struct_grid, line={\"row\": 10})\nxs.plot_grid();","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:41.926453Z","iopub.execute_input":"2024-04-11T01:15:41.927465Z","iopub.status.idle":"2024-04-11T01:15:42.50011Z","shell.execute_reply.started":"2024-04-11T01:15:41.927418Z","shell.execute_reply":"2024-04-11T01:15:42.4988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cuadrícula estructurada con espaciado variable entre filas y columnas.","metadata":{}},{"cell_type":"code","source":"dx = dy = 5000\n# dx = dy = 200\n\n# creando celdas de transición\nmultiplier = 1.175\ntransition = 20000.0\nncells = 7\nsmoothr = [\n    transition * (multiplier - 1.0) / (multiplier ** float(ncells) - 1.0)\n]\nfor i in range(ncells - 1):\n    smoothr.append(smoothr[i] * multiplier)\nsmooth = smoothr.copy()\nsmooth.reverse()","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:42.501594Z","iopub.execute_input":"2024-04-11T01:15:42.50199Z","iopub.status.idle":"2024-04-11T01:15:42.509111Z","shell.execute_reply.started":"2024-04-11T01:15:42.501956Z","shell.execute_reply":"2024-04-11T01:15:42.507914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"construir la cuadrícula estructurada con espaciado variable entre filas y columnas","metadata":{}},{"cell_type":"code","source":"dx = 9 * [5000] + smooth + 15 * [1666.66666667] + smoothr + 14 * [5000]\ndy = 4 * [5000] + smooth + 12 * [1666.66666667] + smoothr + 4 * [5000]\n\nnlay = 1\nncol = len(dx)\nnrow = len(dy)\n\ntop = np.ones((nrow, ncol)) * 1000.0\nbotm = np.ones((nlay, nrow, ncol)) * -100.0\n\ndelr = np.array(dx)\ndelc = np.array(dy)\nstruct_vrc_grid = StructuredGrid(\n    nlay=nlay,\n    delr=delr,\n    delc=delc,\n    xoff=0.0,\n    yoff=0.0,\n    top=top,\n    botm=botm,\n)\nset_idomain(struct_vrc_grid, boundary_polygon)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:42.511042Z","iopub.execute_input":"2024-04-11T01:15:42.511477Z","iopub.status.idle":"2024-04-11T01:15:42.938998Z","shell.execute_reply.started":"2024-04-11T01:15:42.511436Z","shell.execute_reply":"2024-04-11T01:15:42.937586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_sg_vrc = fine_topo.resample_to_grid(\n    struct_vrc_grid,\n    band=fine_topo.bands[0],\n    method=\"linear\",\n    extrapolate_edges=True,\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:42.942115Z","iopub.execute_input":"2024-04-11T01:15:42.943044Z","iopub.status.idle":"2024-04-11T01:15:43.441471Z","shell.execute_reply.started":"2024-04-11T01:15:42.942991Z","shell.execute_reply":"2024-04-11T01:15:43.440102Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ixs = flopy.utils.GridIntersect(struct_vrc_grid, method=\"structured\")\ncellids = []\nfor sg in sgs:\n    v = ixs.intersect(LineString(sg), sort_by_cellid=True)\n    cellids += v[\"cellids\"].tolist()\nintersection_sg_vrc = np.zeros(struct_vrc_grid.shape[1:])\nfor loc in cellids:\n    intersection_sg_vrc[loc] = 1","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:43.44335Z","iopub.execute_input":"2024-04-11T01:15:43.444205Z","iopub.status.idle":"2024-04-11T01:15:43.521724Z","shell.execute_reply.started":"2024-04-11T01:15:43.444154Z","shell.execute_reply":"2024-04-11T01:15:43.520085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=figsize)\nax = fig.add_subplot()\npmv = flopy.plot.PlotMapView(modelgrid=struct_vrc_grid)\nax.set_aspect(\"equal\")\npmv.plot_array(top_sg_vrc)\npmv.plot_array(\n    intersection_sg_vrc,\n    masked_values=[\n        0,\n    ],\n    alpha=0.2,\n    cmap=\"Reds_r\",\n)\ncg = pmv.contour_array(\n    top_sg_vrc, levels=levels, linewidths=0.3, colors=\"0.75\"\n)\npmv.plot_inactive()\n\nax.plot(bp[:, 0], bp[:, 1], \"k-\")\nfor sg in sgs:\n    sa = np.array(sg)\n    ax.plot(sa[:, 0], sa[:, 1], \"b-\", alpha=0.2)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:43.522951Z","iopub.execute_input":"2024-04-11T01:15:43.523275Z","iopub.status.idle":"2024-04-11T01:15:43.973666Z","shell.execute_reply.started":"2024-04-11T01:15:43.523246Z","shell.execute_reply":"2024-04-11T01:15:43.97256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Red de refinamiento de red local","metadata":{}},{"cell_type":"code","source":"from flopy.utils.lgrutil import Lgr\n\nnlayp = 1\ndx = 5000\nnrowp = int(Ly / dx)\nncolp = int(Lx / dx)\ndelrp = dx\ndelcp = dx\n\ntopp = np.ones((nrowp, ncolp)) * 1000.0\nbotmp = np.ones((nlayp, nrowp, ncolp)) * -100.0\n\nidomainp = np.ones((nlayp, nrowp, ncolp), dtype=int)\nidomainp[0, 8:12, 13:18] = 0\n\nncpp = 3\nncppl = [1]\n\nlgr = Lgr(\n    nlayp,\n    nrowp,\n    ncolp,\n    delrp,\n    delcp,\n    topp,\n    botmp,\n    idomainp,\n    ncpp=ncpp,\n    ncppl=ncppl,\n    xllp=0.0,\n    yllp=0.0,\n)\n\ndelr = np.array(ncolp * [dx])\ndelc = np.array(nrowp * [dx])\nstruct_gridp = StructuredGrid(\n    nlay=1, delr=delr, delc=delc, idomain=idomainp, top=topp, botm=botmp\n)\nset_idomain(struct_gridp, boundary_polygon)\n\ndelr_child, delc_child = lgr.get_delr_delc()\nxoff, yoff = lgr.get_lower_left()\nnrowc, ncolc = delc_child.shape[0], delr_child.shape[0]\ntopc = np.ones((nrowc, ncolc)) * 1000.0\nbotmc = np.ones((nlayp, nrowc, ncolc)) * -100.0\nstruct_gridc = StructuredGrid(\n    delr=delr_child,\n    delc=delc_child,\n    xoff=xoff,\n    yoff=yoff,\n    idomain=idomainp,\n    top=topc,\n    botm=botmc,\n)\n\nextent_child = struct_gridc.extent\n\nnested_grid = [struct_gridp, struct_gridc]","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:43.975237Z","iopub.execute_input":"2024-04-11T01:15:43.975988Z","iopub.status.idle":"2024-04-11T01:15:44.158867Z","shell.execute_reply.started":"2024-04-11T01:15:43.975944Z","shell.execute_reply":"2024-04-11T01:15:44.157653Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_ngp = fine_topo.resample_to_grid(\n    struct_gridp,\n    band=fine_topo.bands[0],\n    method=\"linear\",\n    extrapolate_edges=True,\n)\ntop_ngc = fine_topo.resample_to_grid(\n    struct_gridc,\n    band=fine_topo.bands[0],\n    method=\"linear\",\n    extrapolate_edges=True,\n)\ntop_nested_grid = [top_ngp, top_ngc]","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:44.16036Z","iopub.execute_input":"2024-04-11T01:15:44.161253Z","iopub.status.idle":"2024-04-11T01:15:45.145156Z","shell.execute_reply.started":"2024-04-11T01:15:44.161217Z","shell.execute_reply":"2024-04-11T01:15:45.143506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ixs = flopy.utils.GridIntersect(struct_gridp, method=\"structured\")\ncellids = []\nfor sg in sgs:\n    v = ixs.intersect(LineString(sg), sort_by_cellid=True)\n    cellids += v[\"cellids\"].tolist()\nintersection_ngp = np.zeros(struct_gridp.shape[1:])\nfor loc in cellids:\n    intersection_ngp[loc] = 1\nintersection_ngp[idomainp[0] == 0] = 0\n\nixs = flopy.utils.GridIntersect(struct_gridc, method=\"structured\")\ncellids = []\nfor sg in sgs:\n    v = ixs.intersect(LineString(sg), sort_by_cellid=True)\n    cellids += v[\"cellids\"].tolist()\nintersection_ngc = np.zeros(struct_gridc.shape[1:])\nfor loc in cellids:\n    intersection_ngc[loc] = 1\n\nintersection_nested_grid = [intersection_ngp, intersection_ngc]","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:45.146988Z","iopub.execute_input":"2024-04-11T01:15:45.147525Z","iopub.status.idle":"2024-04-11T01:15:45.220674Z","shell.execute_reply.started":"2024-04-11T01:15:45.147481Z","shell.execute_reply":"2024-04-11T01:15:45.219507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=figsize)\nax = fig.add_subplot()\npmv = flopy.plot.PlotMapView(modelgrid=struct_gridp, extent=extent)\npmv.plot_inactive()\npmv.plot_array(top_ngp, vmin=vmin, vmax=vmax)\npmv.plot_array(\n    intersection_nested_grid[0],\n    masked_values=[\n        0,\n    ],\n    alpha=0.2,\n    cmap=\"Reds_r\",\n)\ncgp = pmv.contour_array(top_ngp, levels=levels, linewidths=0.3, colors=\"0.75\")\npmv.plot_inactive(zorder=100)\nax.set_aspect(\"equal\")\n\npmvc = flopy.plot.PlotMapView(modelgrid=struct_gridc, ax=ax, extent=extent)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:45.222625Z","iopub.execute_input":"2024-04-11T01:15:45.223021Z","iopub.status.idle":"2024-04-11T01:15:45.690972Z","shell.execute_reply.started":"2024-04-11T01:15:45.222984Z","shell.execute_reply":"2024-04-11T01:15:45.689495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pmvc.plot_grid()\npmvc.plot_array(top_ngc, vmin=vmin, vmax=vmax)\npmvc.plot_array(\n    intersection_nested_grid[1],\n    masked_values=[\n        0,\n    ],\n    alpha=0.2,\n    cmap=\"Reds_r\",\n)\ncgc = pmvc.contour_array(top_ngc, levels=levels, linewidths=0.3, colors=\"0.75\")\n\nax.plot(bp[:, 0], bp[:, 1], \"k-\")\nfor sg in sgs:\n    sa = np.array(sg)\n    ax.plot(sa[:, 0], sa[:, 1], \"b-\")","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:45.692522Z","iopub.execute_input":"2024-04-11T01:15:45.692907Z","iopub.status.idle":"2024-04-11T01:15:45.726866Z","shell.execute_reply.started":"2024-04-11T01:15:45.692875Z","shell.execute_reply":"2024-04-11T01:15:45.725323Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Cuadrículas no estructuradas\n\ncrear un polígono del área de refinamiento común para la cuadrícula de quadtree","metadata":{}},{"cell_type":"code","source":"lgr_poly = [\n    [\n        (extent_child[0], extent_child[2]),\n        (extent_child[0], extent_child[3]),\n        (extent_child[1], extent_child[3]),\n        (extent_child[1], extent_child[2]),\n        (extent_child[0], extent_child[2]),\n    ]\n]","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:45.728957Z","iopub.execute_input":"2024-04-11T01:15:45.729356Z","iopub.status.idle":"2024-04-11T01:15:45.736457Z","shell.execute_reply.started":"2024-04-11T01:15:45.729319Z","shell.execute_reply":"2024-04-11T01:15:45.735018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Grillado de Quadtree ","metadata":{}},{"cell_type":"code","source":"sim = flopy.mf6.MFSimulation()\ngwf = gwf = flopy.mf6.ModflowGwf(sim)\ndx = dy = 5000.0\nnr = int(Ly / dy)\nnc = int(Lx / dx)\ndis6 = flopy.mf6.ModflowGwfdis(\n    gwf,\n    nrow=nr,\n    ncol=nc,\n    delr=dy,\n    delc=dx,\n)\ng = Gridgen(gwf.modelgrid, model_ws=temp_path)\nadpoly = [\n    [[(1000, 1000), (3000, 1000), (3000, 2000), (1000, 2000), (1000, 1000)]]\n]\nadpoly = boundary_polygon + [boundary_polygon[0]]\nadpoly = [[adpoly]]\ng.add_refinement_features([lgr_poly], \"polygon\", 2, range(1))\n\nrefine_line = sgs\ng.add_refinement_features(refine_line, \"line\", 2, range(1))\ng.build(verbose=False)\n\ngridprops_vg = g.get_gridprops_vertexgrid()\nquadtree_grid = flopy.discretization.VertexGrid(**gridprops_vg)\nset_idomain(quadtree_grid, boundary_polygon)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:45.738104Z","iopub.execute_input":"2024-04-11T01:15:45.738518Z","iopub.status.idle":"2024-04-11T01:15:47.689584Z","shell.execute_reply.started":"2024-04-11T01:15:45.738483Z","shell.execute_reply":"2024-04-11T01:15:47.688446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_qg = fine_topo.resample_to_grid(\n    quadtree_grid,\n    band=fine_topo.bands[0],\n    method=\"linear\",\n    extrapolate_edges=True,\n)\n\nixs = flopy.utils.GridIntersect(quadtree_grid, method=\"vertex\")\ncellids = []\nfor sg in sgs:\n    v = ixs.intersect(LineString(sg), sort_by_cellid=True)\n    cellids += v[\"cellids\"].tolist()\nintersection_qg = np.zeros(quadtree_grid.shape[1:])\nfor loc in cellids:\n    intersection_qg[loc] = 1","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:47.691055Z","iopub.execute_input":"2024-04-11T01:15:47.691504Z","iopub.status.idle":"2024-04-11T01:15:48.593809Z","shell.execute_reply.started":"2024-04-11T01:15:47.691471Z","shell.execute_reply":"2024-04-11T01:15:48.592373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=figsize)\nax = fig.add_subplot()\npmv = flopy.plot.PlotMapView(modelgrid=quadtree_grid)\npmv.plot_array(top_qg, ec=\"0.75\")\npmv.plot_array(\n    intersection_qg,\n    masked_values=[\n        0,\n    ],\n    alpha=0.2,\n    cmap=\"Reds_r\",\n)\ncg = pmv.contour_array(top_qg, levels=levels, linewidths=0.3, colors=\"white\")\npmv.plot_inactive(zorder=100)\nax.set_aspect(\"equal\")\n\nax.plot(bp[:, 0], bp[:, 1], \"k-\")\nfor sg in sgs:\n    sa = np.array(sg)\n    ax.plot(sa[:, 0], sa[:, 1], \"b-\")","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:48.595688Z","iopub.execute_input":"2024-04-11T01:15:48.597119Z","iopub.status.idle":"2024-04-11T01:15:49.178824Z","shell.execute_reply.started":"2024-04-11T01:15:48.597069Z","shell.execute_reply":"2024-04-11T01:15:49.177603Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Rejilla triangular","metadata":{}},{"cell_type":"code","source":"nodes = []\nfor sg in sgs:\n    sg_densify = densify_geometry(sg, 2000)\n    nodes.extend(sg_densify)\nfor x in struct_gridc.get_xvertices_for_layer(0)[0, :]:\n    for y in struct_gridc.get_yvertices_for_layer(0)[:, 0]:\n        nodes.append((x, y))\nnodes = np.array(nodes)\n","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:49.180704Z","iopub.execute_input":"2024-04-11T01:15:49.181121Z","iopub.status.idle":"2024-04-11T01:15:49.2162Z","shell.execute_reply.started":"2024-04-11T01:15:49.181084Z","shell.execute_reply":"2024-04-11T01:15:49.214917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tri = Triangle(maximum_area=5000 * 5000, angle=30, nodes=nodes, model_ws=temp_path)\n# tri = Triangle(maximum_area=1000 * 1000, angle=30, nodes=nodes, model_ws=temp_path)\npoly = bp\ntri.add_polygon(poly)\ntri.build(verbose=False)\n\n\ncell2d = tri.get_cell2d()\nvertices = tri.get_vertices()\n\ntop = np.ones(tri.ncpl) * 1000.0\nbotm = np.ones((1, tri.ncpl)) * -100.0\nidomain = np.ones((1, tri.ncpl), dtype=int)\n\ntriangular_grid = VertexGrid(\n    vertices=vertices,\n    cell2d=cell2d,\n    idomain=idomain,\n    nlay=1,\n    ncpl=tri.ncpl,\n    top=top,\n    botm=botm,\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:49.218857Z","iopub.execute_input":"2024-04-11T01:15:49.220238Z","iopub.status.idle":"2024-04-11T01:15:49.319432Z","shell.execute_reply.started":"2024-04-11T01:15:49.220187Z","shell.execute_reply":"2024-04-11T01:15:49.318142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_tg = fine_topo.resample_to_grid(\n    triangular_grid,\n    band=fine_topo.bands[0],\n    method=\"linear\",\n    extrapolate_edges=True,\n)\n\nixs = flopy.utils.GridIntersect(triangular_grid)  # , method=\"vertex\")\ncellids = []\nfor sg in sgs:\n    v = ixs.intersect(\n        LineString(sg),\n        return_all_intersections=True,\n        keepzerolengths=False,\n        sort_by_cellid=True,\n    )\n    cellids += v[\"cellids\"].tolist()\nintersection_tg = np.zeros(triangular_grid.shape[1:])\nfor loc in cellids:\n    intersection_tg[loc] = 1","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:49.320856Z","iopub.execute_input":"2024-04-11T01:15:49.321199Z","iopub.status.idle":"2024-04-11T01:15:50.065097Z","shell.execute_reply.started":"2024-04-11T01:15:49.321169Z","shell.execute_reply":"2024-04-11T01:15:50.063642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=figsize)\nax = fig.add_subplot()\nax.set_aspect(\"equal\")\n\npmv = flopy.plot.PlotMapView(modelgrid=triangular_grid)\n\npmv.plot_array(top_tg, ec=\"0.75\")\npmv.plot_array(\n    intersection_tg,\n    masked_values=[\n        0,\n    ],\n    alpha=0.2,\n    cmap=\"Reds_r\",\n)\ncg = pmv.contour_array(top_tg, levels=levels, linewidths=0.3, colors=\"white\")\nax.clabel(cg, cg.levels, inline=True, fmt=\"%1.0f\", fontsize=10)\n\npmv.plot_inactive(zorder=100)\n\nif True:\n    ax.plot(bp[:, 0], bp[:, 1], \"k-\")\n    for sg in sgs:\n        sa = np.array(sg)\n        ax.plot(sa[:, 0], sa[:, 1], \"b-\")","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:50.066651Z","iopub.execute_input":"2024-04-11T01:15:50.067048Z","iopub.status.idle":"2024-04-11T01:15:50.870964Z","shell.execute_reply.started":"2024-04-11T01:15:50.067014Z","shell.execute_reply":"2024-04-11T01:15:50.869296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Voronoi Grid from the Triangle object\nCrear un objeto vor y VertexGrid a partir del objeto triangular tri","metadata":{}},{"cell_type":"code","source":"vor = VoronoiGrid(tri)\ngridprops = vor.get_gridprops_vertexgrid()\nidomain = np.ones((1, vor.ncpl), dtype=int)\nvoronoi_grid = VertexGrid(**gridprops, nlay=1, idomain=idomain)\n\ntop_vg = fine_topo.resample_to_grid(\n    voronoi_grid,\n    band=fine_topo.bands[0],\n    method=\"linear\",\n    extrapolate_edges=True,\n)\n\nixs = flopy.utils.GridIntersect(voronoi_grid, method=\"vertex\")\ncellids = []\nfor sg in sgs:\n    v = ixs.intersect(\n        LineString(sg),\n        return_all_intersections=True,\n        keepzerolengths=False,\n        sort_by_cellid=True,\n    )\n    cellids += v[\"cellids\"].tolist()\nintersection_vg = np.zeros(voronoi_grid.shape[1:])\nfor loc in cellids:\n    intersection_vg[loc] = 1","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:50.872593Z","iopub.execute_input":"2024-04-11T01:15:50.873007Z","iopub.status.idle":"2024-04-11T01:15:51.78187Z","shell.execute_reply.started":"2024-04-11T01:15:50.87297Z","shell.execute_reply":"2024-04-11T01:15:51.780487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=figsize)\nax = fig.add_subplot()\npmv = flopy.plot.PlotMapView(modelgrid=voronoi_grid)\nax.set_aspect(\"equal\")\npmv.plot_array(top_vg)\npmv.plot_array(\n    intersection_vg,\n    masked_values=[\n        0,\n    ],\n    alpha=0.2,\n    cmap=\"Reds_r\",\n)\npmv.plot_inactive()\nax.plot(bp[:, 0], bp[:, 1], \"k-\")\nfor sg in sgs:\n    sa = np.array(sg)\n    ax.plot(sa[:, 0], sa[:, 1], \"b-\")\n\ncg = pmv.contour_array(top_vg, levels=levels, linewidths=0.3, colors=\"0.75\")","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:51.783724Z","iopub.execute_input":"2024-04-11T01:15:51.784238Z","iopub.status.idle":"2024-04-11T01:15:52.332626Z","shell.execute_reply.started":"2024-04-11T01:15:51.78419Z","shell.execute_reply":"2024-04-11T01:15:52.328977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Ploteo de las 6 cuadriculas\nCrear una polilínea para el área de Refinamiento de cuadrícula local","metadata":{}},{"cell_type":"code","source":"lgr_poly_array = np.array(lgr_poly).squeeze()","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:52.334187Z","iopub.execute_input":"2024-04-11T01:15:52.335399Z","iopub.status.idle":"2024-04-11T01:15:52.341345Z","shell.execute_reply.started":"2024-04-11T01:15:52.335348Z","shell.execute_reply":"2024-04-11T01:15:52.34015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"figwidth = 17.15 / 2.54\nfigheight = 3.66 * (Ly / Lx) * 8.25 / 2.54\ngrids = [\n    struct_grid,\n    struct_vrc_grid,\n    nested_grid,\n    quadtree_grid,\n    triangular_grid,\n    voronoi_grid,\n    None,\n]\ntopo_grids = [top_sg, top_sg_vrc, top_nested_grid, top_qg, top_tg, top_vg]\nextent = (0, 180000, 0, 100000)\ncbar_axis = [0.0, 0.1, 0.25, 0.9]\n\nwith styles.USGSMap():\n    fig = plt.figure(figsize=(figwidth, figheight), constrained_layout=True)\n    gs = gridspec.GridSpec(ncols=16, nrows=17, figure=fig)\n    axs = [fig.add_subplot(gs[:5, :8])]\n    axs.append(fig.add_subplot(gs[:5, 8:]))\n    axs.append(fig.add_subplot(gs[5:10, :8]))\n    axs.append(fig.add_subplot(gs[5:10, 8:]))\n    axs.append(fig.add_subplot(gs[10:15, :8]))\n    axs.append(fig.add_subplot(gs[10:15, 8:]))\n    axs.append(fig.add_subplot(gs[15:, :]))\n\n    for idx, ax in enumerate(axs[:-1]):\n        g = grids[idx]\n        t = topo_grids[idx]\n        if g is not None:\n            if isinstance(g, list):\n                g = g[0]\n                t = t[0]\n\n            pmv = flopy.plot.PlotMapView(modelgrid=g, ax=ax)\n\n            v = pmv.plot_array(t)\n            pmv.plot_grid(**grid_dict)\n            cg = pmv.contour_array(t, **contour_dict)\n            pmv.ax.clabel(cg, cg.levels, **clabel_dict)\n            pmv.plot_inactive(color_noflow=\"gray\", zorder=100)\n\n            if isinstance(grids[idx], list):\n                gg = grids[idx]\n                tt = topo_grids[idx]\n                for g, t in zip(gg[1:], tt[1:]):\n                    pmvc = flopy.plot.PlotMapView(\n                        modelgrid=g, ax=ax, extent=extent\n                    )\n                    pmvc.plot_array(top_ngc, vmin=vmin, vmax=vmax)\n                    pmvc.plot_grid(**grid_dict)\n                    cgc = pmvc.contour_array(top_ngc, **contour_dict)\n\n            # plot lgr polyline\n            ax.plot(\n                lgr_poly_array[:, 0],\n                lgr_poly_array[:, 1],\n                zorder=101,\n                **refinement_dict,\n            )\n\n            ax.set_aspect(\"equal\")\n            ax.set_axisbelow(False)\n\n            ax.set_xlim(extent[0], extent[1])\n            ax.set_xticks(np.arange(0, 200000, 50000))\n            if idx in (4, 5):\n                ax.set_xticklabels(np.arange(0, 200, 50))\n                ax.set_xlabel(\"x position (km)\")\n            else:\n                ax.set_xticklabels([])\n\n            ax.set_ylim(extent[2], extent[3])\n            ax.set_yticks(np.arange(0, 150000, 50000))\n            if idx in (0, 2, 4):\n                ax.set_yticklabels(np.arange(0, 150, 50))\n                ax.set_ylabel(\"y position (km)\")\n            else:\n                ax.set_yticklabels([])\n\n            styles.heading(ax, idx=idx)\n            if True:\n                ax.plot(bp[:, 0], bp[:, 1], \"k-\", lw=2.0)\n                for sg in sgs:\n                    sa = np.array(sg)\n                    ax.plot(sa[:, 0], sa[:, 1], **river_dict)\n\n    # legend\n    ax = axs[-1]\n    xy0 = (-100, -100)\n    ax.set_xlim(0, 1)\n    ax.set_ylim(0, 1)\n    ax.set_axis_off()\n\n    ax.axhline(\n        xy0[0],\n        color=\"black\",\n        lw=2,\n        label=\"Basin boundary\",\n    )\n    ax.axhline(\n        xy0[0],\n        **river_dict,\n        label=\"River\",\n    )\n    ax.axhline(\n        xy0[0],\n        color=contour_color,\n        lw=0.5,\n        ls=\"--\",\n        label=\"Elevation contour\",\n    )\n    ax.axhline(\n        xy0[0],\n        label=\"Grid refinement area\",\n        **refinement_dict,\n    )\n    ax.axhline(\n        xy0[0],\n        marker=\"s\",\n        mec=\"0.5\",\n        mfc=\"white\",\n        markeredgewidth=0.25,\n        lw=0,\n        label=\"Grid cell\",\n    )\n    ax.axhline(\n        xy0[0],\n        marker=\"s\",\n        mec=\"0.5\",\n        mfc=\"gray\",\n        markeredgewidth=0.25,\n        lw=0,\n        label=\"Inactive area\",\n    )\n    styles.graph_legend(\n        ax,\n        ncol=3,\n        loc=\"lower center\",\n        labelspacing=0.5,\n        columnspacing=0.6,\n        handletextpad=0.3,\n    )\n\n    # color de barra\n    ax = fig.add_subplot(gs[15:, 14:])\n    ax.set_xlim(0, 1)\n    ax.set_ylim(0, 1)\n    ax.set_axis_off()\n    cax = ax.inset_axes(\n        cbar_axis,\n    )\n    # cax.set_axisbelow(False)\n    cbar = plt.colorbar(\n        v,\n        orientation=\"vertical\",\n        cax=cax,\n        ticks=[25, 50, 75, 100],\n    )\n    cbar.ax.tick_params(\n        labelsize=5,\n        labelcolor=\"black\",\n        color=\"black\",\n        length=9,\n        pad=2,\n    )\n    cbar.ax.set_title(\n        \"Elevation (m)\", pad=2.5, loc=\"center\", fontdict=font_dict\n    )\n    \n    ","metadata":{"execution":{"iopub.status.busy":"2024-04-11T01:15:52.343464Z","iopub.execute_input":"2024-04-11T01:15:52.344151Z","iopub.status.idle":"2024-04-11T01:15:56.097295Z","shell.execute_reply.started":"2024-04-11T01:15:52.344106Z","shell.execute_reply":"2024-04-11T01:15:56.096067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Ploteo del rio intersectado para las 6 cuadriculas","metadata":{}},{"cell_type":"code","source":"figwidth = 17.15 / 2.54\nfigheight = 3.66 * (Ly / Lx) * 8.25 / 2.54\ngrids = [\n    struct_grid,\n    struct_vrc_grid,\n    nested_grid,\n    quadtree_grid,\n    triangular_grid,\n    voronoi_grid,\n    None,\n]\nintersections = [\n    intersection_sg,\n    intersection_sg_vrc,\n    intersection_nested_grid,\n    intersection_qg,\n    intersection_tg,\n    intersection_vg,\n    None,\n]\nextent = list(extent_child)\nde = 10000.0\nextent[0] -= 2.0 * de\nextent[1] += 2.5 * de\nextent[2] -= de\nextent[3] += de\n\n\nwith styles.USGSMap():\n    fig = plt.figure(figsize=(figwidth, figheight), constrained_layout=True)\n    gs = gridspec.GridSpec(ncols=2, nrows=16, figure=fig)\n    axs = [fig.add_subplot(gs[:5, 0])]\n    axs.append(fig.add_subplot(gs[:5, 1]))\n    axs.append(fig.add_subplot(gs[5:10, 0]))\n    axs.append(fig.add_subplot(gs[5:10, 1]))\n    axs.append(fig.add_subplot(gs[10:15, 0]))\n    axs.append(fig.add_subplot(gs[10:15, 1]))\n    axs.append(fig.add_subplot(gs[15:, :]))\n\n    for idx, ax in enumerate(axs[:-1]):\n        g = grids[idx]\n        t = intersections[idx]\n        if g is not None:\n            if isinstance(g, list):\n                g = g[0]\n                t = t[0]\n\n            pmv = flopy.plot.PlotMapView(modelgrid=g, ax=ax, extent=extent)\n\n            v = pmv.plot_array(t, masked_values=(0,), cmap=intersection_cmap)\n            pmv.plot_grid(**grid_dict)\n\n            pmv.plot_inactive(color_noflow=\"gray\", zorder=100)\n\n            if isinstance(grids[idx], list):\n                gg = grids[idx]\n                tt = intersections[idx]\n                for g, t in zip(gg[1:], tt[1:]):\n                    pmvc = flopy.plot.PlotMapView(\n                        modelgrid=g, ax=ax, extent=extent\n                    )\n                    pmvc.plot_array(\n                        t, masked_values=(0,), cmap=intersection_cmap\n                    )\n                    pmvc.plot_grid(**grid_dict)\n\n            # plot lgr polyline\n            ax.plot(\n                lgr_poly_array[:, 0],\n                lgr_poly_array[:, 1],\n                zorder=101,\n                **refinement_dict,\n            )\n\n            ax.set_aspect(\"equal\")\n            ax.set_axisbelow(False)\n\n            ax.set_xlim(extent[0], extent[1])\n            ax.set_xticks(np.arange(50000, 120000, 10000))\n            if idx in (4, 5):\n                ax.set_xticklabels(np.arange(50, 120, 10))\n                ax.set_xlabel(\"x position (km)\")\n            else:\n                ax.set_xticklabels([])\n\n            ax.set_ylim(extent[2], extent[3])\n            ax.set_yticks(np.arange(35000, 70000, 10000))\n            if idx in (0, 2, 4):\n                ax.set_yticklabels(np.arange(35, 75, 10))\n                ax.set_ylabel(\"y position (km)\")\n            else:\n                ax.set_yticklabels([])\n\n            styles.heading(ax, idx=idx)\n            if True:\n                ax.plot(bp[:, 0], bp[:, 1], \"k-\", lw=2.0)\n                for sg in sgs:\n                    sa = np.array(sg)\n                    ax.plot(sa[:, 0], sa[:, 1], **river_dict)\n\n    # legend\n    ax = axs[-1]\n    xy0 = (-100, -100)\n    ax.set_xlim(0, 1)\n    ax.set_ylim(0, 1)\n    ax.set_axis_off()\n\n    ax.axhline(xy0[0], **river_dict, label=\"River\")\n    ax.axhline(\n        xy0[0],\n        label=\"Grid refinement area\",\n        **refinement_dict,\n    )\n    ax.axhline(\n        xy0[0],\n        marker=\"s\",\n        mec=\"0.5\",\n        mfc=\"white\",\n        markeredgewidth=0.25,\n        lw=0,\n        label=\"Grid cell\",\n    )\n    ax.axhline(\n        xy0[0],\n        marker=\"s\",\n        mfc=intersection_rgba,\n        mec=\"0.5\",\n        markeredgewidth=0.25,\n        lw=0,\n        label=\"Intersected cell\",\n    )\n    styles.graph_legend(\n        ax,\n        ncol=4,\n        loc=\"lower center\",\n        labelspacing=0.5,\n        columnspacing=0.6,\n        handletextpad=0.3,\n    )","metadata":{"tags":[],"execution":{"iopub.status.busy":"2024-04-11T01:15:56.099229Z","iopub.execute_input":"2024-04-11T01:15:56.099615Z","iopub.status.idle":"2024-04-11T01:15:58.598975Z","shell.execute_reply.started":"2024-04-11T01:15:56.099581Z","shell.execute_reply":"2024-04-11T01:15:58.598059Z"},"trusted":true},"execution_count":null,"outputs":[]}]}