{ "cells": [ { "cell_type": "markdown", "id": "5af6d7ab-3c43-4800-a66a-cbd36d5b9a4e", "metadata": {}, "source": [ "## Example 2" ] }, { "cell_type": "code", "execution_count": 2, "id": "treated-cabinet", "metadata": {}, "outputs": [ { "ename": "ModuleNotFoundError", "evalue": "No module named 'cartopy'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", "Input \u001b[0;32mIn [2]\u001b[0m, in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyplot\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mplt\u001b[39;00m\n\u001b[1;32m 6\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtri\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mtri\u001b[39;00m\n\u001b[0;32m----> 7\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mcartopy\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mcrs\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mccrs\u001b[39;00m\n\u001b[1;32m 9\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mnp\u001b[39;00m\n", "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'cartopy'" ] } ], "source": [ "from math import ceil\n", "import pystare as ps\n", "\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "import matplotlib.tri as tri\n", "import cartopy.crs as ccrs\n", "\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "id": "indonesian-duncan", "metadata": {}, "outputs": [], "source": [ "def shiftarg_lon(lon):\n", " \"If lon is outside +/-180, then correct back.\"\n", " if(lon>180):\n", " return ((lon + 180.0) % 360.0)-180.0\n", " else:\n", " return lon\n", "\n", "def triangulate1(lats,lons):\n", " \"Prepare data for tri.Triangulate.\"\n", " print('triangulating1...')\n", " intmat=[]\n", " npts=int(len(lats)/3)\n", " k=0\n", " for i in range(npts):\n", " intmat.append([k,k+1,k+2])\n", " k=k+3\n", " for i in range(len(lons)):\n", " lons[i] = shiftarg_lon(lons[i])\n", " print('triangulating1 done.') \n", " return lons,lats,intmat\n", "\n", "def plot1(lon,lat,lons,lats,triang,c0='r',c1='b',transf=None,lw=1):\n", " if(lon is not None):\n", " x=np.zeros([lon.size+1],dtype=np.double);x[:-1]=lon[:];x[-1]=lon[0]\n", " y=np.zeros([lat.size+1],dtype=np.double); y[:-1]=lat[:]; y[-1]=lat[0]\n", " ax.plot(x,y,True,transform=transf,c=c0)\n", " plt.triplot(triang,c1+'-',transform=transf,lw=lw,markersize=3)\n", " plt.scatter(lons,lats,s=10,c=c1,transform=ccrs.PlateCarree())\n", " return\n", "\n", "def make_hull(lat0,lon0,resolution0):\n", " hull0 = ps.to_hull_range_from_latlon(lat0,lon0,resolution0)\n", " lath0,lonh0,lathc0,lonhc0 = ps.to_vertices_latlon(hull0)\n", " lons0,lats0,intmat0 = ps.triangulate(lath0,lonh0)\n", " triang0 = tri.Triangulation(lons0,lats0,intmat0)\n", " return lats0,lons0,triang0,hull0" ] }, { "cell_type": "code", "execution_count": null, "id": "coordinated-debut", "metadata": {}, "outputs": [], "source": [ "resolution = 5\n", "# resolution = 4\n", "\n", "resolution0 = resolution\n", "lat0 = np.array([ 10, 5, 60,70], dtype=np.double)\n", "lon0 = np.array([-30,-20,60,10], dtype=np.double)\n", "lats0,lons0,triang0,hull0 = make_hull(lat0,lon0,resolution0)\n", "print('hull0: ',len(hull0))\n", "\n", "resolution1 = resolution\n", "lat1 = np.array([10, 20, 30, 20 ], dtype=np.double)\n", "lon1 = np.array([-60, 60, 60, -60], dtype=np.double)\n", "lats1,lons1,triang1,hull1 = make_hull(lat1,lon1,resolution1)\n", "print('hull1: ',len(hull1))\n", "\n", "if True:\n", " intersected = np.full([1000],-1,dtype=np.int64)\n", " # intersected = ps.intersect(hull0,hull1,multiresolution=False)\n", " intersected = ps.intersection(hull0, hull1, multi_resolution=True)\n", " # intersected = ps.intersect(hull0,hull1,multiresolution=True)\n", " # print('hull0: ',[hex(i) for i in hull0])\n", " # print('hull1: ',[hex(i) for i in hull1])\n", " # ps._intersect_multiresolution(hull0,hull1,intersected)\n", " # print('intersected: ',len(intersected))\n", " # print('np.min: ',np.amin(intersected))\n", " # print('intersected: ',[hex(i) for i in intersected])\n", " # The following are for _intersect_multiresolution's results\n", " # endarg = np.argmax(intersected < 0)\n", " # intersected = intersected[:endarg]\n", " # intersected = ps.intersect(hull0,hull1)\n", " print('intersected: ',len(intersected))\n", " lati,loni,latci,lonci = ps.to_vertices_latlon(intersected)\n", " lonsi,latsi,intmati = ps.triangulate(lati,loni)\n", " triangi = tri.Triangulation(lonsi,latsi,intmati)\n", "\n", "\n", "# Set up the projection and transformation\n", "proj = ccrs.PlateCarree()\n", "# proj = ccrs.Robinson()\n", "# proj = ccrs.Geodesic()\n", "# proj = ccrs.Mollweide()\n", "transf = ccrs.Geodetic()\n", "# transf = ccrs.PlateCarree()\n", "plt.figure()\n", "plt.subplot(projection=proj,transform=transf)\n", "ax = plt.axes(projection=proj,transform=transf)\n", "ax.set_global()\n", "ax.coastlines()\n", "print('graphics-0')\n", "plot1(lon0,lat0,lons0,lats0,triang0,c0='r',c1='b',transf=transf)\n", "print('graphics-1')\n", "plot1(lon1,lat1,lons1,lats1,triang1,c0='g',c1='c',transf=transf)\n", "print('graphics-2')\n", "plot1(None,None,lonsi,latsi,triangi,c0='y',c1='r',transf=transf,lw=4)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "stone-moses", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.4" } }, "nbformat": 4, "nbformat_minor": 5 }