diff --git a/README.md b/README.md index de24794..e16e52e 100644 --- a/README.md +++ b/README.md @@ -3,4 +3,4 @@ Available on [http://voila-gallery.pascalbugnion.net](http://voila-gallery.pascalbugnion.net/voila/render/index). -To play with the code, open this in [Binder](https://mybinder.org/v2/gh/pbugnion/voila-gallery/master?urlpath=voila%2Ftree). +To play with the code, open this in [Binder](https://mybinder.org/v2/gh/erdc/voila-gallery/master?urlpath=voila%2Ftree). diff --git a/index.ipynb b/index.ipynb index 52c6520..ab5ec84 100644 --- a/index.ipynb +++ b/index.ipynb @@ -14,7 +14,9 @@ "\n", "- [country-indicators](country-indicators/index) shows how to integrate Matplotlib plots using output widgets.\n", "- [gaussian-density](gaussian-density/index) shows how to build a dashboard around a bqplot plot.\n", - "- [stl-display](stl-display/index) shows how to read STL files and display them with [ipyvolume](https://github.com/maartenbreddels/ipyvolume)." + "- [stl-display](stl-display/index) shows how to read STL files and display them with [ipyvolume](https://github.com/maartenbreddels/ipyvolume).\n", + "- [mesh-generation-ipymesh](mesh-generation-ipymesh/index) shows how to construct a Planar Straight Line Graph with [ipymesh](https://github.com/erdc/ipymesh) and triangulate it with [triangle](https://www.cs.cmu.edu/~quake/triangle.html).\n", + "- [mesh-generation-ipyleaflet](mesh-generation-ipyleaflet/index) shows how to construct a Planar Straight Line Graph with [ipyleaflet](https://github.com/jupyter-widgets/ipymesh) and triangulate it with [triangle](https://www.cs.cmu.edu/~quake/triangle.html), adaptively, to approximate complex bathymetry/terrain from the GEBCO global grid.\n" ] }, { @@ -41,7 +43,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.4" + "version": "3.7.2" } }, "nbformat": 4, diff --git a/mesh-generation-ipyleaflet/index.ipynb b/mesh-generation-ipyleaflet/index.ipynb new file mode 100644 index 0000000..eafeb62 --- /dev/null +++ b/mesh-generation-ipyleaflet/index.ipynb @@ -0,0 +1,256 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this app we:\n", + "* Create a very simple [Planar Straight Line Graph](https://en.wikipedia.org/wiki/Planar_straight-line_graph) with the [ipyleaflet](https://github.com/jupygter-widgets/ipyleaflet) widget. In fact, it's just a polygon, and the rest of the app can only handle a single polygon at a time right now--so trash your existing polygon before creating another one. \n", + "* Pull the terrain/bathymetry data for the area of interest from a global database. The notebook has code to use [GEBCO 2019](https://www.gebco.net). In fact, the app is pulling from a small grid around the initial coordinates to keep the docker container small.\n", + "* Generate a triangular mesh with the [triangle](https://www.cs.cmu.edu/~quake/triangle.html) mesh generator using the [triangle python bindings](https://rufat.be/triangle/index.html) and adaptively refine it to resolve the terrain data.\n", + "* Visualize the resulting mesh in several ways.\n", + "Find the code [here](https://github.com/pbugnion/voila-gallery/blob/master/mesh-generation-ipyleaflet/index.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib widget\n", + "import ipywidgets\n", + "from ipyleaflet import (\n", + " Map,\n", + " Polygon,\n", + " DrawControl,\n", + " WMSLayer\n", + ")\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import triangle as tr\n", + "from matplotlib import tri as mpltr\n", + "from scipy import interpolate as scipy_interpolate\n", + "import ipyvolume as ipv\n", + "import ipympl\n", + "\n", + "#initial map center and zoom level\n", + "center = [34.6252978589571, -77.34580993652344]\n", + "vertices=np.array([[-77.521791, 34.582935],\n", + " [-77.380628, 34.356501],\n", + " [-77.11704 , 34.561348],\n", + " [-77.321914, 34.735944]])\n", + "zoom = 9\n", + "#use the GEBCO tile layers in the map widget--this is not necessary, just for viz\n", + "wms = WMSLayer(url='http://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv', \n", + " layers='GEBCO_LATEST', \n", + " )\n", + "#create a drawing control \n", + "dc = DrawControl(polyline={},#\n", + " marker={},#{'shapeOptions': {'color': '#0000FF'}},\n", + " rectangle={},#{'shapeOptions': {'color': '#0000FF'}},\n", + " circle={},#{'shapeOptions': {'color': '#0000FF'}},\n", + " circlemarker={}#,\n", + " )\n", + "\n", + "m = Map(center=center, zoom=zoom, layers=(wms,), controls=(dc,))" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "plt.ioff()\n", + "gridfig=plt.figure(num=0)\n", + "gridfig.set_label(' ')\n", + "#import h5py\n", + "#f = h5py.File('GEBCO_2019.nc','r')\n", + "f = np.load(\"terrain_slice_save.npz\")\n", + "def plotGrid(button=None):\n", + " global vertices,dc, j_min, j_max, i_min, i_max, lonr, latr, gridfig\n", + " #get the vertices of the polygon\n", + " if dc.last_draw['geometry'] != None:\n", + " vertices=np.array(dc.last_draw['geometry']['coordinates'][0][:-1])\n", + " lonr = [vertices[:,0].min(),vertices[:,0].max()]\n", + " latr = [vertices[:,1].min(),vertices[:,1].max()]\n", + " j_min = np.abs(f['lon'][:] - lonr[0]).argmin() - 1\n", + " j_max = np.abs(f['lon'][:] - lonr[1]).argmin() + 1\n", + " i_min = np.abs(f['lat'][:] - latr[0]).argmin() - 1\n", + " i_max = np.abs(f['lat'][:] - latr[1]).argmin() + 1\n", + " plt.figure(gridfig.number)\n", + " plt.ion()\n", + " plt.clf()\n", + " plt.contourf(f['lon'][j_min:j_max],\n", + " f['lat'][i_min:i_max],\n", + " f['elevation'][i_min:i_max,\n", + " j_min:j_max])\n", + " plt.colorbar()\n", + " plt.title('')\n", + "plotGrid()\n", + "terrain_button = ipywidgets.Button(\n", + " description='Get Terrain',\n", + " disabled=False,\n", + " button_style='', # 'success', 'info', 'warning', 'danger' or ''\n", + " tooltip='Pull the terrain for the area of interest'\n", + ")\n", + "terrain_button.on_click(plotGrid)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "plt.ioff()\n", + "meshfig=plt.figure(num=1,frameon=True)\n", + "meshfig.set_label(' ')\n", + "surffig=ipv.figure()\n", + "def genMesh(button=None):\n", + " global vertices,meshfig\n", + " nv = vertices.shape[0]\n", + " vertexFlags=[1 for i in range(nv)]\n", + " segments=[]\n", + " segmentFlags=[]\n", + " for i in range(nv):\n", + " segments.append([i,(i+1)%nv])\n", + " segmentFlags.append(1)\n", + " nx = j_max - j_min\n", + " ny = i_max - i_min\n", + " bathy = np.zeros((nx*ny,3),'d')\n", + " for i in range(ny):\n", + " bathy[i*nx:(i+1)*nx,0] = f['lon'][j_min:j_max]\n", + " bathy[i*nx:(i+1)*nx,1] = f['lat'][i_min+i]\n", + " bathy[:,2] = np.reshape(f['elevation'][i_min:i_max,j_min:j_max], nx*ny)\n", + "\n", + " he = min((lonr[1]-lonr[0])/40.,(latr[1]-latr[0])/40.0)\n", + " bathyGridDim = (ny,nx)\n", + " x = bathy[:bathyGridDim[1],0]\n", + " y = bathy[:bathyGridDim[0]*bathyGridDim[1]:bathyGridDim[1],1]\n", + " z = bathy[:,2].reshape(bathyGridDim).transpose()\n", + " bathyInterpolant = scipy_interpolate.RectBivariateSpline(x,y,z,kx=1,ky=1)\n", + " pslg = dict(vertices=vertices,\n", + " vertex_attributes=vertexFlags,\n", + " segments=segments,\n", + " segment_markers=segmentFlags)\n", + " mesh = tr.triangulate(pslg,'pq30Da{0:f}'.format(0.5*he**2))\n", + " area_max = 0.5*he**2\n", + " mplmesh=mpltr.Triangulation(mesh['vertices'][:,0],\n", + " mesh['vertices'][:,1],\n", + " mesh['triangles'])\n", + " not_converged=True\n", + " tol=0.5\n", + " while not_converged:\n", + " mplmesh=mpltr.Triangulation(mesh['vertices'][:,0],\n", + " mesh['vertices'][:,1],\n", + " mesh['triangles'])\n", + " trifinder = mpltr.TrapezoidMapTriFinder(mplmesh)\n", + " zNodes = bathyInterpolant.ev(mesh['vertices'][:,0],\n", + " mesh['vertices'][:,1])\n", + " interpolator = mpltr.LinearTriInterpolator(mplmesh, zNodes, trifinder=trifinder)\n", + " zInterp = interpolator(bathy[:,0],bathy[:,1])\n", + " xE = bathy[np.abs(bathy[:,2] - zInterp) > tol,0]\n", + " yE = bathy[np.abs(bathy[:,2] - zInterp) > tol,1]\n", + " e = trifinder(xE,yE)\n", + " area=np.ones((mesh['triangles'].shape[0],),'d')*0.5*he**2\n", + " area_max *= 0.5\n", + " area[e[e>=0]] = area_max\n", + " mesh['triangle_max_area']=area\n", + " mesh2 = tr.triangulate(mesh,'rpq30Da'.format(0.5*he**2))\n", + " err = np.abs(bathy[:,2] - zInterp) > tol\n", + " #print(\"L_infty_error\",np.abs(bathy[:,2] - zInterp).max())\n", + " not_converged = err.any()\n", + " mesh=mesh2\n", + " fig=plt.figure(meshfig.number)\n", + " plt.ion()\n", + " plt.clf()\n", + " plt.triplot(mplmesh)\n", + " fig = ipv.figure(surffig)\n", + " if len(fig.meshes):\n", + " del fig.meshes[-1]\n", + " ipv.style.use('minimal')\n", + " ipv.plot_trisurf(x=mesh['vertices'][:,0],\n", + " y=mesh['vertices'][:,1],\n", + " z=zNodes,\n", + " triangles=mesh['triangles'],\n", + " color='grey')\n", + " #,color=color, wireframe=wireframe)\n", + " ipv.xlim(lonr[0], lonr[1])\n", + " ipv.ylim(latr[0], latr[1])\n", + " ipv.zlim(zNodes.min()*10,zNodes.max()*10)\n", + " mesh3d = ipv.gcc()\n", + "button = ipywidgets.Button(\n", + " description='Generate Mesh',\n", + " disabled=False,\n", + " button_style='', # 'success', 'info', 'warning', 'danger' or ''\n", + " tooltip='Regenerate the mesh'\n", + ")\n", + "genMesh()\n", + "button.on_click(genMesh)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "4a9030355a504927b16682c93ce8b80b", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Tab(children=(Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'max_zoom': 19, 'attri…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from ipywidgets import Tab,VBox\n", + "app=Tab(children=[m,\n", + " VBox(children=[terrain_button, gridfig.canvas]),\n", + " VBox(children=[button, meshfig.canvas]),\n", + " surffig])\n", + "app.set_title(0,\"Map\")\n", + "app.set_title(1,\"Terrain Data\")\n", + "app.set_title(2,\"Mesh\")\n", + "app.set_title(3,\"Terrain Surface\")\n", + "display(app)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "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.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/mesh-generation-ipyleaflet/terrain_slice_save.npz b/mesh-generation-ipyleaflet/terrain_slice_save.npz new file mode 100644 index 0000000..807a53b Binary files /dev/null and b/mesh-generation-ipyleaflet/terrain_slice_save.npz differ diff --git a/mesh-generation-ipymesh/index.ipynb b/mesh-generation-ipymesh/index.ipynb new file mode 100644 index 0000000..d004133 --- /dev/null +++ b/mesh-generation-ipymesh/index.ipynb @@ -0,0 +1,173 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this app we:\n", + "* Edit a [Planar Straight Line Graph](https://en.wikipedia.org/wiki/Planar_straight-line_graph) with the [erdc/ipymesh](https://github.com/erdc/ipymesh) widget\n", + "* Generate a triangular mesh with the [triangle](https://www.cs.cmu.edu/~quake/triangle.html) mesh generator using the [triangle python bindings](https://rufat.be/triangle/index.html)\n", + "\n", + "Find the code [here](https://github.com/pbugnion/voila-gallery/blob/master/mesh-generation/index.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from ipymesh import PSLGEditor\n", + "editor = PSLGEditor(vertices=[(0.0,0.0),\n", + " (0.0,1.0),\n", + " (1.0,1.0),\n", + " (1.0,0.0)],\n", + " vertexFlags=[1,1,2,2],\n", + " segments=[(0,1),\n", + " (1,2),\n", + " (2,3),\n", + " (3,0)],\n", + " segmentFlags=[1,1,2,1],\n", + " regions=[(0.5,0.5)],\n", + " regionFlags=[1],\n", + " regionTypes=[1,2,0],\n", + " boundaryTypes=[1,2,0],\n", + " Lx=1.05,\n", + " Ly=1.05,\n", + " x0=-0.15,\n", + " y0=-0.15,\n", + " width=500,\n", + " height=500)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import ipywidgets\n", + "import triangle as tr\n", + "import matplotlib.pyplot as plt\n", + "from bqplot import (\n", + " LinearScale, OrdinalColorScale, Lines, Axis, ColorAxis, Figure\n", + ")\n", + "fig=Figure()\n", + "def genMesh(b=None):\n", + " global editor,fig\n", + " pslg = dict(vertices=editor.vertices,\n", + " vertex_attributes=editor.vertexFlags,\n", + " segments=editor.segments,\n", + " segment_markers=editor.segmentFlags)\n", + " if editor.holes:\n", + " pslg['holes'] = editor.holes\n", + " if editor.regions:\n", + " pslg['regions'] = [[r[0],r[1],rF, 0.1] for (r,rF) in \n", + " zip(editor.regions,editor.regionFlags)]\n", + " mesh = tr.triangulate(pslg,'pq')\n", + " size = 20\n", + " x_data = [[x1,x2] for x1,x2 in zip(mesh['vertices'][mesh['segments'][:,0],0],\n", + " mesh['vertices'][mesh['segments'][:,1],0])]\n", + " y_data = [[y1,y2] for y1,y2 in zip(mesh['vertices'][mesh['segments'][:,0],1],\n", + " mesh['vertices'][mesh['segments'][:,1],1])]\n", + " colorsList = [ci[0] for ci in mesh['segment_markers']]\n", + " edges=set()\n", + " for t in mesh['triangles']:\n", + " for n in range(3):\n", + " e = [t[n],t[(n+1)%3]]\n", + " e.sort()\n", + " edges.add(tuple(e))\n", + " for e in edges:\n", + " x_data.append([mesh['vertices'][e[0],0],mesh['vertices'][e[1],0]])\n", + " y_data.append([mesh['vertices'][e[0],1],mesh['vertices'][e[1],1]])\n", + " colorsList.append(0)\n", + " assert(len(colorsList) == len(x_data))\n", + " x_sc = LinearScale()\n", + " y_sc = LinearScale()\n", + " c_sc = OrdinalColorScale(domain=[0,1,2])\n", + " line = Lines(x=x_data, \n", + " y=y_data, \n", + " color=colorsList,\n", + " scales={'x': x_sc, \n", + " 'y': y_sc, \n", + " 'color': c_sc},\n", + " stroke_width=3, \n", + " display_legend=False)\n", + " ax_x = Axis(scale=x_sc, grid_lines='solid', label='x')\n", + " ax_y = Axis(scale=y_sc, orientation='vertical', tick_format='0.2f',\n", + " grid_lines='solid', label='y')\n", + " ax_c = ColorAxis(scale=c_sc, side='right', visible=True)\n", + "\n", + " fig.marks=[line]\n", + " fig.axes=[ax_x, ax_y, ax_c]\n", + " fig.max_aspect_ratio = 1\n", + " fig.min_aspect_ratio = 1\n", + "genMesh()\n", + "button = ipywidgets.Button(\n", + " description='Generate Mesh',\n", + " disabled=False,\n", + " button_style='', # 'success', 'info', 'warning', 'danger' or ''\n", + " tooltip='Regenerate the mesh'\n", + ")\n", + "button.on_click(genMesh)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "71dc594058894e88923a15802c126eee", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Tab(children=(VBox(children=(Button(description='Generate Mesh', style=ButtonStyle(), tooltip='Regenerate the …" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from ipywidgets import Tab,VBox\n", + "app=Tab(children=[VBox(children=[button, editor]),\n", + " fig])\n", + "app.set_title(0,\"Graph Editor\")\n", + "app.set_title(1,\"Mesh\")\n", + "display(app)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "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.7.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/requirements.txt b/requirements.txt index 9b4db98..50eda36 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ matplotlib==2.2.3 -ipywidgets +ipywidgets==7.4.2 +ipympl==0.2.1 voila==0.0.13 pandas==0.23.4 bqplot==0.11.5 @@ -7,3 +8,6 @@ scipy==1.1.0 ipyvolume==0.5.1 numpy-stl==2.10.1 ipyupload==0.1.3 +ipymesh==0.1.0 +triangle==20190115.2 +ipyleaflet==0.10.2