{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "cvt_florida.ipynb\n", "\n", "Discussion: Compute a CVT using the shape of the state of Florida.\n", "\n", "Licensing: This code is distributed under the GNU LGPL license.\n", " \n", "Modified: 10 November 2016\n", "\n", "Author: John Burkardt, Lukas Bystricky" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using matplotlib backend: agg\n" ] } ], "source": [ "# Import necessary libraries and set plot option\n", "%matplotlib\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import math\n", "import scipy.spatial as spatial" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# A CVT of Florida #\n", "\n", "We can approximate the shape of Florida using a polygon, whose\n", "vertices are given as longitude and latitude values. Then we\n", "can pick the locations of a few Florida cities as initial generators and\n", "apply our CVT iteration. \n", "\n", "We will look at:\n", "* how to read a file containing polygon information, and display it.\n", "* another way to determine if a point is in a polygon.\n", "* another way to uniformly randomly sample a polygon.\n", "* how to use the longitude and latitude of some Florida cities as starting generators.\n", "* how to compute one step of the CVT;\n", "* how to run the CVT and display the final configuration." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Longitude and Latitudes for Florida Polygon Vertices#\n", "\n", "The file \"florida_shape.txt\" is available from the web page \n", "http://people.sc.fsu.edu/~jburkardt/classes/urop_2016/urop_2016.html\n", "\n", "The Florida polygon consists of 570 vertices. The data file lists the\n", "longitude and latitude of each vertex, which you can treat as though\n", "these were X and Y coordinates. The file repeats the first line as the\n", "last line (#571) which is a plus if you are drawing the polygon.\n", "\n", "We need to read this file, and make a plot of the points to convince ourselves\n", "that we read it correctly." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0 -81.5000 30.7200\n", " 1 -81.4800 30.6800\n", " 2 -81.4500 30.7000\n", " 3 -81.4300 30.6800\n", " 4 -81.4500 30.5700\n", " 5 -81.4300 30.5200\n", " 6 -81.4700 30.5500\n", " 7 -81.4800 30.6200\n", " 8 -81.5000 30.6200\n", " 9 -81.5000 30.5500\n", "...\n", " 560 -81.7800 30.7700\n", " 561 -81.7300 30.7700\n", " 562 -81.7200 30.7500\n", " 563 -81.6800 30.7300\n", " 564 -81.6500 30.7300\n", " 565 -81.6300 30.7300\n", " 566 -81.6200 30.7200\n", " 567 -81.5500 30.7000\n", " 568 -81.5200 30.7200\n", " 569 -81.5000 30.7200\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Code to read and display longitude and latitude for Florida polygon vertices\n", "\n", "def florida_shape_read ( ):\n", "\n", " import numpy as np\n", "\n", " input = open ( 'florida_shape.txt', 'r' )\n", "\n", " lon = np.zeros ( 570 )\n", " lat = np.zeros ( 570 )\n", "\n", " i = 0\n", " for line in input:\n", "\n", " if ( line[0] == '#' ):\n", " continue\n", " else:\n", " data = line.split ( )\n", " lon[i] = data[0]\n", " lat[i] = data[1]\n", " i = i + 1\n", "\n", " input.close ( )\n", "\n", " return lon, lat\n", "#\n", "# Test the function.\n", "#\n", "lon, lat = florida_shape_read ( )\n", "for i in range ( 0, 10 ):\n", " print ( ' %3d %8.4f %8.4f' % ( i, lon[i], lat[i] ) ) \n", "print ( '...')\n", "for i in range ( 560, 570 ):\n", " print ( ' %3d %8.4f %8.4f' % ( i, lon[i], lat[i] ) ) \n", "\n", "plt.plot ( lon, lat, 'b-' )\n", "plt.axis ( 'Equal')\n", "plt.title ( 'The outline of Florida')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Polygon contains point? #\n", "\n", "In the polygon notebook, we showed that we used triangulation\n", "to determine if a polygon contained a point. When dealing with\n", "Florida, we will prefer another method, which is made available to\n", "you in a Python file:\n", " \n", "http://people.sc.fsu.edu/~jburkardt/classes/urop_2016/polygon_contains_point.py\n", " \n", "You should copy this and put it in the same directory as your notebook,\n", "and invoke it using the command\n", "\n", " from polygon_contains_point import polygon_contains_point\n", " \n", "The function has the form\n", "```python\n", "def polygon_contains_point ( x, y, px, py ):\n", " ***\n", " return inside\n", "```\n", "where the point is (x,y), the coordinates of the vertices of the \n", "polygon are in the arrays px and py, and the result \"inside\" is \n", "True if the polygon contains the point.\n", "\n", "Test this function using the florida longitude and latitude data\n", "for the polygon coordinates px and py, and the following points:\n", " \n", " Point Lon Lat Where?\n", " 0 -87 26 (In the ocean!)\n", " 1 -81 30 (In the ocean!)\n", " 2 -84 30.5 (Near Tallahassee\n", " 3 -81 26 (Near Miami)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0 -87.0000 26.0000 False\n", " 1 -81.0000 30.0000 False\n", " 2 -84.0000 30.5000 True\n", " 3 -81.0000 26.0000 True\n" ] } ], "source": [ "# Test code for polygon contains point #\n", "\n", "from polygon_contains_point import polygon_contains_point\n", "\n", "x = np.array ( [ -87.0, -81.0, -84.0, -81.0 ] )\n", "y = np.array ( [ 26.0, 30.0, 30.5, 26.0 ] )\n", "\n", "for i in range ( 0, 4 ):\n", " inside = polygon_contains_point ( x[i], y[i], lon, lat )\n", " print ( ' %2d %8.4f %8.4f %s' % ( i, x[i], y[i], inside ) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Acceptance/Rejection Sampling of the Florida Polygon #\n", "\n", "In the polygon notebook, we showed that we could use triangulation\n", "to sample the polygon uniformly. But because Florida is a complicated\n", "polygon, we will try a simpler method. \n", "\n", "To get a uniformly random point from Florida:\n", "* Surround Florida by a rectangular box;\n", "* Get a random point in that box (easy);\n", "* If the Florida polygon does not contain that point, go back one step;\n", "\n", "A box that contains Florida has the lower left corner (-88,25) and\n", "the upper right corner (-80,31).\n", "\n", "To get a sample point, then, we get unit random numbers r1 and r2 and \n", "compute:\n", "\n", " x = -88 + 8 * r1\n", " y = 25 + 6 * r2" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Code for acceptance/rejection sampling of the Florida Polygon\n", "\n", "def florida_sample ( m ):\n", " xy = np.zeros ( [ m, 2 ] )\n", " i = 0\n", " while ( i < m ):\n", " r1 = np.random.rand ( )\n", " r2 = np.random.rand ( )\n", " x = -88.0 + 8.0 * r1\n", " y = 25.0 + 6.0 * r2\n", " inside = polygon_contains_point ( x, y, lon, lat )\n", " if ( inside ):\n", " xy[i,0] = x\n", " xy[i,1] = y\n", " i = i + 1\n", " return xy\n", "#\n", "# Test the code.\n", "#\n", "m = 100\n", "xy = florida_sample ( m )\n", "plt.plot ( lon, lat, 'b-' )\n", "plt.plot ( xy[:,0], xy[:,1], 'ro' )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# CVT step for Florida #\n", "\n", "Now we can write the CVT step for Florida.\n", "\n", "Given N current generators G, and using M sample points, we want to\n", "estimate the N centroids C of the corresponding Voronoi regions.\n", "\n", "Use the following seven Florida cities for G:\n", "\n", " -87.20, 30.43 Pensacola\n", " -84.25, 30.46 Tallahassee\n", " -81.66, 30.34 Jacksonville\n", " -82.48, 27.97 Tampa\n", " -82.32, 29.65 Gainesville\n", " -81.30, 28.42 Orlando\n", " -80.21, 25.78 Miami" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[-86.61426552 30.71168878]\n", " [-84.54031681 30.32754013]\n", " [-81.70440116 30.21489764]\n", " [-82.12693031 27.66841477]\n", " [-82.50193766 29.61688588]\n", " [-81.17685053 28.12548008]\n", " [-80.83812526 26.16262678]]\n" ] } ], "source": [ "# cvt_step_florida()\n", "#\n", "def cvt_step_florida ( g, m ):\n", " import numpy as np\n", " n = g.shape[0]\n", " s = florida_sample ( m )\n", " ni = np.zeros ( n )\n", " c = np.zeros ( [ n, 2 ])\n", " for i in range ( 0, m ):\n", " k = -1\n", " d = np.Inf\n", " for j in range ( 0, n ):\n", " dj = np.linalg.norm ( s[i,:] - g[j,:] )\n", " if ( dj < d ):\n", " d = dj\n", " k = j\n", " ni[k] = ni[k] + 1\n", " c[k,:] = c[k,:] + s[i,:]\n", " \n", " for i in range ( 0, n ):\n", " c[i,:] = c[i,:] / float ( ni[i] )\n", " \n", " return c\n", "#\n", "# Test the function using 7 Florida cities to start with.\n", "#\n", "import numpy as np\n", "g = np.array ( [ \\\n", " [ -87.20, 30.43 ], \\\n", " [ -84.25, 30.46 ], \\\n", " [ -81.66, 30.34 ], \\\n", " [ -82.48, 27.97 ], \\\n", " [ -82.32, 29.65 ], \\\n", " [ -81.30, 28.42 ], \\\n", " [ -80.21, 25.78 ] ] )\n", "\n", "m = 1000\n", "c = cvt_step_florida ( g, m )\n", "print ( c )" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# CVT Iteration for Florida #\n", "\n", "Once we have a CVT step function, the CVT iteration is easy. Try taking 10 steps,\n", "and plot the locations of G before the iteration, and then at the end.\n", "\n", "You may also be able to plot the Voronoi diagram on top of the plot of the cities." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Code for CVT Iteration for Florida\n", "#\n", "def cvt_florida ( n, m, step_num ):\n", "\n", " g = florida_sample ( n )\n", "\n", "# vor = spatial.Voronoi ( g )\n", "# spatial.voronoi_plot_2d ( vor )\n", " plt.plot ( g[:,0], g[:,1], 'go' )\n", " plt.plot ( lon, lat, 'b-' )\n", " plt.axis ( 'Equal' )\n", " plt.title ( 'Before CVT iteration')\n", " plt.show ( )\n", "\n", " for step in range ( 0, step_num ):\n", " c = cvt_step_florida ( g, m )\n", " g = c.copy ( )\n", " return g\n", "#\n", "# Test the code.\n", "#\n", "n = 7\n", "m = 5000\n", "step_num = 10\n", "g = cvt_florida ( n, m, step_num )\n", "\n", "#vor = spatial.Voronoi ( g )\n", "#spatial.voronoi_plot_2d ( vor )\n", "plt.plot ( g[:,0], g[:,1], 'ro' )\n", "plt.plot ( lon, lat, 'b-' )\n", "plt.axis ( 'Equal' )\n", "plt.title ( 'After CVT iteration')\n", "plt.show ( )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.8" } }, "nbformat": 4, "nbformat_minor": 0 }