transition_nhc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "86b07cea-2d91-4a4c-ad49-f2233d110fc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ocsmesh\n",
    "import geopandas as gpd\n",
    "from pyproj import CRS,Transformer\n",
    "from shapely import intersection\n",
    "import pandas as pd\n",
    "from stormevents.nhc import VortexTrack\n",
    "from copy import deepcopy\n",
    "from ocsmesh import Hfun, Mesh, JigsawDriver,Geom\n",
    "import numpy as np\n",
    "from scipy.interpolate import griddata\n",
    "from shapely.ops import polygonize, unary_union, transform\n",
    "from shapely.geometry import ( # type: ignore[import]\n",
    "        Polygon, MultiPolygon,\n",
    "        box, GeometryCollection, Point, MultiPoint,\n",
    "        LineString, LinearRing)\n",
    "import utm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3da73ff8-b46b-421f-b022-e847eda3899a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _calculate_mesh_size_function(\n",
    "        buffer_domain,\n",
    "        hires_mesh_clip,\n",
    "        lores_mesh_clip,\n",
    "        buffer_crs\n",
    "    ):\n",
    "\n",
    "    assert buffer_crs == hires_mesh_clip.crs == lores_mesh_clip.crs\n",
    "\n",
    "    # HARDCODED FOR NOW\n",
    "    approx_elem_per_width = 3\n",
    "\n",
    "    msht_hi = deepcopy(hires_mesh_clip)\n",
    "    msht_lo = deepcopy(lores_mesh_clip)\n",
    "\n",
    "    crs = buffer_crs\n",
    "    assert(not buffer_crs.is_geographic)\n",
    "\n",
    "    # calculate mesh size for clipped bits\n",
    "    hfun_hi = Hfun(Mesh(msht_hi))\n",
    "    hfun_hi.size_from_mesh()\n",
    "\n",
    "    hfun_lo = Hfun(Mesh(msht_lo))\n",
    "    hfun_lo.size_from_mesh()\n",
    "\n",
    "    msht_cdt = ocsmesh.utils.triangulate_polygon(\n",
    "        buffer_domain, None, opts='p'\n",
    "    )\n",
    "    msht_cdt.crs = crs\n",
    "\n",
    "    hfun_cdt = Hfun(Mesh(msht_cdt))\n",
    "    hfun_cdt.size_from_mesh()\n",
    "\n",
    "    hfun_cdt_sz = deepcopy(hfun_cdt.msh_t().value) / approx_elem_per_width\n",
    "    msht_cdt.value[:] = hfun_cdt_sz\n",
    "    hfun_rep = Hfun(Mesh(msht_cdt))\n",
    "\n",
    "    mesh_domain_rep = JigsawDriver(\n",
    "        geom=Geom(buffer_domain, crs=crs),\n",
    "        hfun=hfun_rep,\n",
    "        initial_mesh=False\n",
    "    ).run(sieve=0)\n",
    "\n",
    "    msht_domain_rep = deepcopy(mesh_domain_rep.msh_t)\n",
    "#        utils.reproject(msht_domain_rep, crs)\n",
    "\n",
    "    pts_2mesh = np.vstack(\n",
    "        (hfun_hi.msh_t().vert2['coord'], hfun_lo.msh_t().vert2['coord'])\n",
    "    )\n",
    "    val_2mesh = np.vstack(\n",
    "        (hfun_hi.msh_t().value, hfun_lo.msh_t().value)\n",
    "    )\n",
    "    domain_sz_1 = griddata(\n",
    "        pts_2mesh, val_2mesh, msht_domain_rep.vert2['coord'], method='linear'\n",
    "    )\n",
    "    domain_sz_2 = griddata(\n",
    "        pts_2mesh, val_2mesh, msht_domain_rep.vert2['coord'], method='nearest'\n",
    "    )\n",
    "    domain_sz = domain_sz_1.copy()\n",
    "    domain_sz[np.isnan(domain_sz_1)] = domain_sz_2[np.isnan(domain_sz_1)]\n",
    "\n",
    "    msht_domain_rep.value[:] = domain_sz\n",
    "    hfun_interp = Mesh(msht_domain_rep)\n",
    "\n",
    "    return hfun_interp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bb97ea21-74ed-4317-926e-7785ffd3261a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_shape(storm, year):\n",
    "\n",
    "    track = VortexTrack.from_storm_name(\n",
    "        storm, year, file_deck='b', advisories=['BEST']\n",
    "    )\n",
    "\n",
    "    ws34 = track.wind_swaths(wind_speed=34)\n",
    "    # Get single best track\n",
    "    best_ws34 = list(ws34['BEST'].values())[0]\n",
    "\n",
    "    return best_ws34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d138231e-31db-4f45-912d-d2149415dd72",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Felicio.Cassalho\\.conda\\envs\\model_vis\\lib\\site-packages\\stormevents\\nhc\\track.py:1109: RuntimeWarning: invalid value encountered in divide\n",
      "  speeds = pandas.Series(distances / abs(intervals), index=indices)\n"
     ]
    }
   ],
   "source": [
    "shape = get_shape('beryl', 2024)\n",
    "\n",
    "highres = ocsmesh.Mesh.open(r'C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/STOFS-3D-Atlantic_v2.1.gr3', crs=4326)\n",
    "lowres = ocsmesh.Mesh.open(r'C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/WNAT_1km.14', crs=4326)\n",
    "\n",
    "highres_clip = ocsmesh.utils.clip_mesh_by_shape(\n",
    "    highres.msh_t,\n",
    "    shape=shape,\n",
    "    inverse=False,\n",
    "    fit_inside=False,\n",
    "    check_cross_edges=False,\n",
    "    adjacent_layers=0\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "17a6ee99-ba97-4c36-b774-54fd2a819f9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(lowres.msh_t).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/lowres.2dm\", format='2dm', overwrite=True)\n",
    "ocsmesh.Mesh(highres.msh_t).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/highres.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b9972d42-683a-438c-99de-dfa70fd3722e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(highres_clip).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/highres_clip.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e129e492-2c57-4dc6-8f05-26335c24722d",
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_gdf = gpd.GeoDataFrame(\n",
    "    geometry = [shape],crs=4326)\n",
    "# shape_gdf.to_file(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/shape.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8e2a8ac5-6d44-4544-aa88-993c19cad0ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged = ocsmesh.utils.merge_overlapping_meshes(\n",
    "    [lowres.msh_t, highres_clip], adjacent_layers=10\n",
    ")\n",
    "\n",
    "ocsmesh.utils.cleanup_duplicates(merged)\n",
    "ocsmesh.utils.put_id_tags(merged)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "87ae56d2-8acf-4840-bf01-5b8aa7a7eee7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(merged).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/merged.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ff76a9f6-d543-4309-b250-c9f2a81210e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "### create_mesh_from_mesh_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "id": "9fd43e01-1bc0-448b-ab18-8ff4be87c430",
   "metadata": {},
   "outputs": [],
   "source": [
    "carved_mesh = ocsmesh.Mesh.open(r'C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\Pearl_River\\temp/carved_mesh.2dm', crs=4326).msh_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "id": "79fc3adf-7f16-4efd-9078-133c4e0b635f",
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_start=[highres_clip,lowres.msh_t]\n",
    "mesh_1 = carved_mesh\n",
    "mesh_2 = highres_clip\n",
    "hfun_mesh = None\n",
    "crs=CRS.from_epsg(4326)\n",
    "min_int_ang=30\n",
    "buffer_domain = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1aa75509-e15b-4a8c-958a-273be54fc587",
   "metadata": {},
   "outputs": [],
   "source": [
    "domain = gpd.GeoDataFrame(geometry=[ocsmesh.utils.get_mesh_polygons(i) for i in domain_start],crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "dfe3fb15-39ca-4e60-91c3-b4368a1772d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_buffer = gpd.GeoDataFrame(geometry=[i[-1].geometry.buffer(buffer_domain) for i in domain.iterrows()],crs=crs)\n",
    "domain_buffer = domain_buffer.dissolve().explode(index_parts=True)\n",
    "domain_buffer.crs = domain_buffer.estimate_utm_crs()\n",
    "domain_buffer =domain_buffer.loc[domain_buffer['geometry'].area == max(domain_buffer['geometry'].area)]\n",
    "domain_buffer.crs = crs\n",
    "domain_buffer = gpd.GeoDataFrame(geometry=[domain_buffer.union_all()],crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0f008471-a744-4fcd-b971-7b3c07d1369b",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_buffer = domain_buffer.union_all().difference(\n",
    "    gpd.GeoDataFrame(\n",
    "        geometry=[\n",
    "            ocsmesh.utils.get_mesh_polygons(mesh_1),\n",
    "            ocsmesh.utils.get_mesh_polygons(mesh_2),\n",
    "        ],\n",
    "        crs = crs\n",
    "    ).union_all()\n",
    ")\n",
    "gdf_full_buffer = gpd.GeoDataFrame(\n",
    "    geometry = [poly_buffer],crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "07f2cb88-c4e2-42ea-8d97-c924c427877f",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipper = domain.union_all().difference(ocsmesh.utils.get_mesh_polygons(mesh_1)).difference(ocsmesh.utils.get_mesh_polygons(mesh_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b6e12a1-e8c5-4063-b953-f9d96fa031f3",
   "metadata": {},
   "source": [
    "### This is the new function, for creating hfun, this is done outside create_mesh_from_mesh_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8030c55-d46f-4ca5-83e9-3991e0fb2113",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh_1_dc = deepcopy(mesh_1)\n",
    "mesh_2_dc = deepcopymesh_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "6b529d4b-cfa7-42c6-8220-a804fff02b23",
   "metadata": {},
   "outputs": [],
   "source": [
    "utm = ocsmesh.utils.estimate_bounds_utm(clipper.bounds, 4326)\n",
    "# Transform all inputs to UTM:\n",
    "t1 = Transformer.from_crs(4326, utm, always_xy=True)\n",
    "clipper_utm = transform(t1.transform, clipper)\n",
    "ocsmesh.utils.reproject(mesh_2_dc, utm)\n",
    "ocsmesh.utils.reproject(mesh_1_dc, utm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "96dc911d-d890-4dda-9a0b-8cb7b7e615a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun_buffer = _calculate_mesh_size_function(\n",
    "        clipper_utm,\n",
    "        mesh_2_dc,\n",
    "        mesh_1_dc,\n",
    "        mesh_1_dc.crs\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "a6353486-863c-4a76-8a5a-85a68451a0cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(hfun_buffer.msh_t).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/hfun_buffer_notutm.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "cc75b757-4301-4f76-aac1-d543a74e678e",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun=Hfun(hfun_buffer)\n",
    "mesh_buf_apprx = JigsawDriver(\n",
    "            geom=Geom(ocsmesh.utils.get_mesh_polygons(hfun_buffer.msh_t), crs=mesh_2_dc.crs),\n",
    "            hfun=hfun,\n",
    "            initial_mesh=False\n",
    "        ).run(sieve=0)\n",
    "\n",
    "ocsmesh.utils.reproject(mesh_buf_apprx.msh_t, lowres.msh_t.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "a3a04462-899e-43e7-83fa-a01c2c23da97",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(mesh_buf_apprx.msh_t).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/mesh_buf_apprx_domain_final.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4544ce4b-88e9-41ed-836d-3624015a46ad",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3153a48d-647e-4eec-9e7a-c57b19301530",
   "metadata": {},
   "source": [
    "#### end of hfun func"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8051065e-12e2-4be1-98c6-a93a16c8cf30",
   "metadata": {},
   "source": [
    "##### If the uses passes a mesh to be the hfun, instead of creating, that should be ok. Inside the create_mesh_from_mesh_diff we should clip the hfun (given or created) with the buffer.\n",
    "##### So, if larger, then no problem/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "id": "22ff8854-b04c-449d-b1dd-50f2ed4808b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun_mesh = mesh_buf_apprx.msh_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "id": "02402fa3-4f26-4a4a-b267-cbd0df57c869",
   "metadata": {},
   "outputs": [],
   "source": [
    "# hfun_clip = ocsmesh.utils.clip_mesh_by_shape(hfun_mesh, clipper, check_cross_edges=True)\n",
    "hfun_clip = hfun_mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "eea8ff10-cc6e-4068-939e-9784266a6c7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(hfun_clip).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/hfun_clip.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "id": "0c31a21f-b070-463e-9d2c-06d772dde2f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary = np.unique(ocsmesh.utils.get_boundary_edges(hfun_clip))\n",
    "all_nodes = hfun_clip.vert2['coord']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "id": "e3b89b26-cd1b-47d3-aa59-0189610b8edf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# arr = np.array(all_nodes)\n",
    "arr = np.delete(all_nodes, boundary, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "id": "95b4b244-3b82-4dc8-978d-aa2640c6ef0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = gpd.GeoDataFrame(geometry = [MultiPoint(arr)],crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "id": "e702e401-3787-43f5-8c53-e7fe65f8a927",
   "metadata": {},
   "outputs": [],
   "source": [
    "root_dir = r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "589ba150-582f-47e5-9b3b-e553ef9a5a04",
   "metadata": {},
   "outputs": [],
   "source": [
    "points.to_file(root_dir+'points.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "bddd9c12-b27c-46d7-b2a3-2bff1835e702",
   "metadata": {},
   "outputs": [],
   "source": [
    "msht_buffer = ocsmesh.utils.triangulate_polygon(gdf_full_buffer,\n",
    "                                    # min_int_ang=30,\n",
    "                                    aux_pts=arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "id": "7b6aac2b-a6eb-4bb5-926c-31bc37b9dfde",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ocsmesh.Mesh(msht_buffer).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/msht_buffer.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "id": "d1c4091f-9a92-462c-8f3b-1ebcc7b40a6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "msht_buffer = ocsmesh.utils.clip_mesh_by_shape(msht_buffer,clipper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "id": "604247cb-32fa-48e2-ba62-994a5f4e3c9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "msht_buffer.crs = crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "id": "46ec501a-47e2-4177-ab78-d3e62966cb82",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(msht_buffer).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/msht_buffer.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "id": "6a317df5-f9ea-402e-8c7c-32354bc28bac",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged = ocsmesh.utils.merge_neighboring_meshes(msht_buffer,mesh_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "id": "8dc70165-4544-45b0-9798-3911e892bb05",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(merged).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/merged.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d9408c9-c7d3-4ecc-bd34-3ab69032d6fe",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a74bed7-22db-4a36-b70f-bafb6df8234c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b932b151-276a-4f97-97ee-327f350bd36f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95e55408-f586-49d9-a3de-9981c0f66a6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary = np.unique(ocsmesh.utils.get_boundary_edges(mesh_buf_apprx))\n",
    "all_nodes = mesh_buf_apprx.vert2['coord']\n",
    "arr = np.delete(all_nodes, boundary, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7537dfe8-bf36-4f81-8b31-f8f23f558d35",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f62d1141-9154-4cb1-836d-89424cc957d5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f85fc14c-31a0-40ab-93e2-9a326781e5e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipper = clipper.union_all().difference(ocsmesh.utils.get_mesh_polygons(mesh_1)).difference(ocsmesh.utils.get_mesh_polygons(mesh_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22cd3a3e-2a81-4972-870b-bde6264756e6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "223a0aa3-ddf4-4309-ba7b-d4157cbc4827",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d94a690f-ca56-49d1-a513-4ae7d8604562",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b20aa460-c7c6-4a21-ab0d-b09453ca438c",
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_cleanup = gpd.GeoDataFrame(geometry=\n",
    "                    gpd.GeoSeries(intersection(\n",
    "                        ocsmesh.utils.get_mesh_polygons(domain_start[0]),\n",
    "                        ocsmesh.utils.get_mesh_polygons(domain_start[1]),\n",
    "                        ))).dissolve().buffer(buffer_domain*2)\n",
    "domain_cleanup = ocsmesh.utils.remove_holes(domain_cleanup.union_all())\n",
    "domain_cleanup = gpd.GeoDataFrame(geometry=\n",
    "                                 gpd.GeoSeries(domain_cleanup),\n",
    "                                 crs=4326).dissolve()\n",
    "domain = pd.concat([gpd.GeoDataFrame\\\n",
    "                    (geometry=[ocsmesh.utils.get_mesh_polygons(i)\\\n",
    "                               .buffer(buffer_domain,join_style=2)\n",
    "                               ],\\\n",
    "                                crs=crs) for i in domain_start])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d0e2e40-d118-44ad-a5f0-57cd4d935322",
   "metadata": {},
   "outputs": [],
   "source": [
    "domain = domain.dissolve().explode(index_parts=True)\n",
    "domain.crs = domain.estimate_utm_crs()\n",
    "domain =domain.loc[domain['geometry'].area == max(domain['geometry'].area)]\n",
    "domain.crs = crs\n",
    "domain = gpd.GeoDataFrame(geometry=[domain.union_all()],crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff29898c-ec6b-4bcb-b42b-111a626f7226",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_buffer = domain.union_all().difference(\n",
    "    gpd.GeoDataFrame(\n",
    "        geometry=[\n",
    "            ocsmesh.utils.get_mesh_polygons(mesh_1),\n",
    "            ocsmesh.utils.get_mesh_polygons(mesh_2),\n",
    "        ],\n",
    "        crs = crs\n",
    "    ).union_all()\n",
    ")\n",
    "gdf_full_buffer = gpd.GeoDataFrame(\n",
    "    geometry = [poly_buffer],crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9fd3f41b-98c7-4801-b6f9-9587e522ac28",
   "metadata": {},
   "outputs": [],
   "source": [
    "root_dir = r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fad9cd66-dc5b-4041-9dfe-5818a78536b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_full_buffer.to_file(root_dir+\"temp/gdf_full_buffer.shp\") \n",
    "domain.to_file(root_dir+\"temp/domain.shp\") \n",
    "domain_cleanup.to_file(root_dir+\"temp/domain_cleanup.shp\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "d26a364d-cb4a-4815-8a44-30e15c305bab",
   "metadata": {},
   "outputs": [],
   "source": [
    "carved_mesh = ocsmesh.Mesh.open(r'C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\Pearl_River\\temp/carved_mesh.2dm', crs=4326).msh_t\n",
    "highres_clip =ocsmesh.Mesh.open(r'C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/highres_clip.2dm', crs=4326).msh_t\n",
    "lowhres =ocsmesh.Mesh.open(r'C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/lowres.2dm', crs=4326).msh_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "327bab8b-ee3b-4c9c-b15d-a096098256e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_final = gpd.GeoDataFrame(\n",
    "        geometry=[\n",
    "            ocsmesh.utils.get_mesh_polygons(highres_clip),\n",
    "            ocsmesh.utils.get_mesh_polygons(lowhres),\n",
    "        ],\n",
    "        crs = 4326\n",
    "    )\n",
    "domain_final = domain_final.union_all().difference(ocsmesh.utils.get_mesh_polygons(carved_mesh)).difference(ocsmesh.utils.get_mesh_polygons(highres_clip))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7327d610-8475-4515-a1a8-ae6f63c2ec07",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpd.GeoDataFrame(geometry = [domain_final],crs=4326).to_file(root_dir+\"temp/domain_final.shp\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "8a74d642-5fda-4da8-b3ad-9d4f06ef34c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_isotach = domain_final"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "de584db7-08d9-4286-ab71-93596c20f61d",
   "metadata": {},
   "outputs": [],
   "source": [
    "utm = ocsmesh.utils.estimate_bounds_utm(poly_isotach.bounds, 4326)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "c11f9b8f-cff9-4c72-815e-f674e5ffd23e",
   "metadata": {},
   "outputs": [],
   "source": [
    "utm = ocsmesh.utils.estimate_bounds_utm(poly_isotach.bounds, 4326)\n",
    "# Transform all inputs to UTM:\n",
    "t1 = Transformer.from_crs(4326, utm, always_xy=True)\n",
    "poly_isotach = transform(t1.transform, poly_isotach)\n",
    "ocsmesh.utils.reproject(highres_clip, utm)\n",
    "ocsmesh.utils.reproject(carved_mesh, utm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "5f09dafb-7172-43e9-aeea-e3047d8c6503",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun_buffer = _calculate_mesh_size_function(\n",
    "        poly_isotach,\n",
    "        highres_clip,\n",
    "        carved_mesh,\n",
    "        carved_mesh.crs\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "41e28b95-dc3b-4f8b-b5b4-1aa3dc9b1ce5",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun=Hfun(hfun_buffer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "d5350e20-91d3-4bc1-bb3d-858efcddfbc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(hfun_buffer.msh_t).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/hfun_buffer_domain_final.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "f01ecb91-586c-4ef3-a174-6bc88c719261",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun_buffer = ocsmesh.Mesh.open(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/hfun_buffer_domain_final.2dm\", crs=4326).msh_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "795b7670-572e-4259-a44d-1422bcbef59d",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun=Hfun(Mesh(hfun_buffer))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "af2641d0-5819-4bae-bede-316f8b389dc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun=Hfun(Mesh(hfun_buffer))\n",
    "mesh_buf_apprx = JigsawDriver(\n",
    "            geom=Geom(ocsmesh.utils.get_mesh_polygons(hfun_buffer), crs=4326),\n",
    "            hfun=hfun,\n",
    "            initial_mesh=False\n",
    "        ).run(sieve=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "b8110690-e8f1-421e-8cc1-bb6f0eff0717",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(mesh_buf_apprx.msh_t).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/mesh_buf_apprx_domain_final.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "af463e7a-d087-44d0-ad6d-9d7a9135b1ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh_buf_apprx = ocsmesh.Mesh.open(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/mesh_buf_apprx_domain_final.2dm\", crs=4326).msh_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "a1fc20f3-48a1-45b6-9601-fe2a84b81d45",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary = np.unique(ocsmesh.utils.get_boundary_edges(mesh_buf_apprx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "fc724d9d-090e-436e-ae8c-0e2631fe4971",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_nodes = mesh_buf_apprx.vert2['coord']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "91673544-c148-40e7-97ae-8829c3b087ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.delete(all_nodes, boundary, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "55500397-c674-4151-8d3c-52c234b32a5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_buffer_mesh = ocsmesh.utils.triangulate_polygon(shape=gdf_full_buffer,aux_pts=arr,opts='p')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "7f16e6f9-5751-4452-8a4a-5a809b644547",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Felicio.Cassalho\\.conda\\envs\\model_vis\\lib\\site-packages\\ocsmesh\\mesh\\mesh.py:106: UserWarning: Input mesh has no CRS information.\n",
      "  warnings.warn('Input mesh has no CRS information.')\n"
     ]
    }
   ],
   "source": [
    "ocsmesh.Mesh(final_buffer_mesh).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/final_buffer_mesh.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "b7e95701-be18-4e75-bc94-7847fc17df16",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_buffer_mesh_clipped = ocsmesh.utils.clip_mesh_by_shape(final_buffer_mesh,domain_final,adjacent_layers=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "1f0c9291-cdb1-40bd-8189-9174cb798b16",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_buffer_mesh_clipped.crs = mesh_buf_apprx.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "6ea08911-d277-4715-9761-bde50a3d33c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(final_buffer_mesh_clipped).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/final_buffer_mesh_clipped.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "a8f26a8f-687b-47a7-80e3-6ebdee7e84c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_merged = ocsmesh.utils.merge_neighboring_meshes(highres_clip,final_buffer_mesh_clipped,carved_mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "e42f9f5d-96a0-4858-908d-3329c97162a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(final_merged).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/final_merged.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32d43bfe-fa40-4736-833c-5a960dbb8ed9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1c99b2b-d99d-4497-8ce6-d83c03a093e9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "837d0b88-5c88-4217-92dc-b210ed42d8ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_full_buffer = gpd.read_file(root_dir+\"temp/gdf_full_buffer.shp\") \n",
    "domain_cleanup = gpd.read_file(root_dir+\"temp/domain_cleanup.shp\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d83e225f-bbc7-4639-8c84-3c2b9f39c6bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly_isotach = intersection(gdf_full_buffer.union_all(), domain_cleanup.union_all())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4d99ebed-2c86-4573-9b8d-3eac1fb4a11b",
   "metadata": {},
   "outputs": [],
   "source": [
    "utm = ocsmesh.utils.estimate_bounds_utm(poly_isotach.bounds, 4326)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5f808647-db98-403b-9c49-8e087cdd8c77",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform all inputs to UTM:\n",
    "t1 = Transformer.from_crs(4326, utm, always_xy=True)\n",
    "poly_isotach = transform(t1.transform, poly_isotach)\n",
    "ocsmesh.utils.reproject(highres_clip, utm)\n",
    "ocsmesh.utils.reproject(carved_mesh, utm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe8eb914-d744-4996-a726-800b1236fc9b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c21dbde-f794-415d-9ebf-66c8a47165e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun_buffer = _calculate_mesh_size_function(\n",
    "        poly_isotach,\n",
    "        highres_clip,\n",
    "        carved_mesh,\n",
    "        carved_mesh.crs\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "457dc41c-4b30-4879-92ee-875d940e7f98",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c077566-ebcb-4d41-b568-0b1dbfc5bfe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun=Hfun(hfun_buffer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "046afea3-b857-4613-abf2-c9b90adb4d9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun_buffer.msh_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "0038fa80-9b59-4ae3-95cf-0a87d935fc33",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(hfun_buffer.msh_t).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/hfun_buffer.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e1ee7cc2-85cd-440b-9098-5cdf5e805af1",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun_buffer = ocsmesh.Mesh.open(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/hfun_buffer.2dm\", crs=4326).msh_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ceaca66f-b35f-4c9f-a9f5-91caf439a0bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "root_dir = r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/\"\n",
    "gdf_full_buffer = gpd.read_file(root_dir+\"temp/gdf_full_buffer.shp\") \n",
    "domain_cleanup = gpd.read_file(root_dir+\"temp/domain_cleanup.shp\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "c5116ba1-ce00-49a1-ae66-b67f8f6b55b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun_buffer = ocsmesh.utils.clip_mesh_by_shape(hfun_buffer,domain_cleanup.union_all())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "958b27bd-84c6-41cc-a992-c4f8a4c090f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "hfun=Hfun(Mesh(hfun_buffer))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "6ccd2f56-f4ee-4418-8496-fcf97f80ae15",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh_buf_apprx = JigsawDriver(\n",
    "            geom=Geom(ocsmesh.utils.get_mesh_polygons(hfun_buffer), crs=4326),\n",
    "            hfun=hfun,\n",
    "            initial_mesh=False\n",
    "        ).run(sieve=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "86a3fe27-9679-4b17-9d2d-ac7205945001",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(mesh_buf_apprx.msh_t).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/mesh_buf_apprx.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "da43a082-60b0-4d56-b32c-cb5296833083",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary = ocsmesh.utils.get_boundary_edges(mesh_buf_apprx.msh_t).ravel()\n",
    "boundary_removed = ocsmesh.utils.remove_mesh_by_edge(mesh_buf_apprx.msh_t,boundary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "c5dc477d-6845-4e19-af55-2801aa850167",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary = ocsmesh.utils.get_boundary_edges(boundary_removed).ravel()\n",
    "boundary_removed = ocsmesh.utils.remove_mesh_by_edge(boundary_removed,boundary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "f1a3cd7c-56da-4b41-86f7-076a6db8ee28",
   "metadata": {},
   "outputs": [],
   "source": [
    "ocsmesh.Mesh(boundary_removed).write(r\"C:\\Users\\Felicio.Cassalho\\Work\\Python_Development\\Mesh\\end-to-end\\subsetting/temp/boundary_removed.2dm\", format='2dm', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f76d9cb1-2a19-4c56-9616-180e3a8061c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_domain = "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "model_vis",
   "language": "python",
   "name": "model_vis"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}