{"id":2456,"date":"2017-10-12T06:21:54","date_gmt":"2017-10-12T11:21:54","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=2456"},"modified":"2017-10-12T06:21:54","modified_gmt":"2017-10-12T11:21:54","slug":"lorens-excellent-adventure-maps-graphs-and-polygons","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2017\/10\/12\/lorens-excellent-adventure-maps-graphs-and-polygons\/","title":{"rendered":"Loren&#8217;s Excellent Adventure: Maps, Graphs, and Polygons"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p><a href=\"https:\/\/www.mathworks.com\/products\/new_products\/latest_features.html\">R2017b<\/a> was released recently. Today's guest blogger, Matt Tearle, discovered a fun application of one of the <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/release-notes.html\">many new features<\/a> - solving a map-based puzzle.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#d77e1159-fdd4-4427-b9f5-b9a45cd298b9\">A traveling MathWorker problem<\/a><\/li><li><a href=\"#34392944-d393-4873-8151-5e78cb017670\">Enter MATLAB<\/a><\/li><li><a href=\"#984a7405-edd6-491e-a323-12f3149f581a\">Problem: who is my neighbor?<\/a><\/li><li><a href=\"#f61b2c99-b43b-4e09-9674-26d62c20a87c\">Solution: the fat of the land<\/a><\/li><li><a href=\"#8e6d4ebf-e568-4117-a1fb-7be89984d444\">Visualizing the results<\/a><\/li><li><a href=\"#f06aac2c-7f31-4aea-b5b8-f4ce0e01d8cd\">Where to next?<\/a><\/li><\/ul><\/div><h4>A traveling MathWorker problem<a name=\"d77e1159-fdd4-4427-b9f5-b9a45cd298b9\"><\/a><\/h4><p><a href=\"http:\/\/www.cartalk.com\/\">\"Car Talk\"<\/a> recently featured a <a href=\"http:\/\/www.cartalk.com\/puzzler\/bobos-dont-look-back-tour\">brain teaser submitted by a former MathWorker<\/a>. The short, slightly modified version is: Loren has to drive from Delaware to every state in the contiguous US, visiting each state exactly once; Loren's boss - let's call him Jack - offers to meet her in the last state, which turns out to be his birthplace; but how can Jack know what the last state will be, without waiting to hear Loren's itinerary?<\/p><p>I'm about to give away the answer, so if you want to figure it out yourself, do so now. It's OK, I'll wait.<\/p><p>Done? Great. Hopefully you realized that there's one and only one state that borders only one other state. That state has to be the beginning or end of any path. (You can't visit any state more than once so once you enter from the single neighboring state, you're stuck.) The beginning is Delaware (which borders several other states), so this unique \"degree-one\" state must be the end of Loren's epic journey.<\/p><h4>Enter MATLAB<a name=\"34392944-d393-4873-8151-5e78cb017670\"><\/a><\/h4><p>As usual for me, the really interesting aspect of this problem is: can I solve it in MATLAB? This is a classic \"traveling salesman\" problem, so the graph capabilities recently added to MATLAB should be of use. But the first problem is how to build the graph that represents the connections between the states. It turns out that another new feature will help.<\/p><p>A new variable type for representing 2-D polygons was introduced in R2017b. The <tt>polyshape<\/tt> function takes vectors representing <i>x<\/i> and <i>y<\/i> coordinates and returns a polygon with those points as its vertices. It so happens that Mapping Toolbox comes with a data file containing the latitudes and longitudes of the borders of the US states - just what I need. With the <tt>shaperead<\/tt> function I can import this data in the form of a structure array, with fields for the name of the state and the latitude\/longitude of the border:<\/p><pre class=\"codeinput\">states = shaperead(<span class=\"string\">'usastatehi.shp'<\/span>);\r\n<span class=\"comment\">% Remove Alaska &amp; Hawaii, and Washington DC<\/span>\r\nstates([2,11,51]) = [];\r\n<\/pre><p>For a single state, it's straightforward to extract the arrays in the <tt>X<\/tt> (longitude) and <tt>Y<\/tt> (latitude) fields and use these as the inputs to <tt>polyshape<\/tt>:<\/p><pre class=\"codeinput\">x = states(1).X;\r\ny = states(1).Y;\r\np = polyshape(x,y)\r\n<\/pre><pre class=\"codeoutput\">p = \r\n  polyshape with properties:\r\n\r\n      Vertices: [518&times;2 double]\r\n    NumRegions: 2\r\n      NumHoles: 0\r\n<\/pre><p>Now that the state is represented as a polygon, I can use the functions that can work with polygon, such as plotting:<\/p><pre class=\"codeinput\">plot(p)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_01.png\" alt=\"\"> <p>But to really get the benefit of working with polygons, I'd like to get an array of them, where each element is a polygon representing one state. I could extract the <i>x<\/i> and <i>y<\/i> components for each state in a loop, but any time I see myself about to apply an operation to each element of an array, I think \"time for <tt>arrayfun<\/tt>!\".<\/p><pre class=\"codeinput\">struct2poly = @(s) polyshape(s.X,s.Y); <span class=\"comment\">% define my struct -&gt; polyshape operation<\/span>\r\np = arrayfun(struct2poly,states) <span class=\"comment\">% apply to all elements of the struct STATES<\/span>\r\n<\/pre><pre class=\"codeoutput\">p = \r\n  48&times;1 polyshape array with properties:\r\n\r\n    Vertices\r\n    NumRegions\r\n    NumHoles\r\n<\/pre><p>Because this is MATLAB, functions like <tt>plot<\/tt> should work on arrays of polygons, just as well as scalar polygons, right?<\/p><pre class=\"codeinput\">plot(p)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_02.png\" alt=\"\"> <p>Isn't that nice? Each polygon is plotted as its own shaded shape. If we zoom in, we can see the solution to the puzzle.<\/p><pre class=\"codeinput\">axis([-75 -66 42 48])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_03.png\" alt=\"\"> <h4>Problem: who is my neighbor?<a name=\"984a7405-edd6-491e-a323-12f3149f581a\"><\/a><\/h4><p>By eye we can see how states share borders, but how can we get MATLAB to determine this in an automated way? Let's start with the simplest problem: given two states (polygons), can we determine if they share a border or not?<\/p><pre class=\"codeinput\">me = p(17);\r\nma = p(19);\r\nnh = p(27);\r\nsubplot(1,2,1)\r\nplot([ma nh])\r\ntitle(<span class=\"string\">'Shared border'<\/span>)\r\nsubplot(1,2,2)\r\nplot([ma me])\r\ntitle(<span class=\"string\">'No shared border'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_04.png\" alt=\"\"> <p>The polygons have a <tt>Vertices<\/tt> property that contains the border locations. So a simplistic approach would be to see if the two polygons have any vertices in common. However, this has a lot of potential problems. Firstly, two shapes can share a boundary without sharing vertices:<\/p><pre class=\"codeinput\">a = polyshape([0 1 1 0],[0 0 3 3]);\r\nb = polyshape([1 2 2 1],[1 1 2 2]);\r\nsubplot(1,1,1)\r\nplot([a b])\r\nintersect(a.Vertices,b.Vertices,<span class=\"string\">'rows'<\/span>)\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n  0&times;2 empty double matrix\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_05.png\" alt=\"\"> <p>Secondly, even if the survey points all align, there's always the potential for floating-point precision issues: how is 31&deg;41'59\" represented?<\/p><pre class=\"codeinput\">lat1 = 31.699722222222222\r\nlat2 = 31 + 41\/60 + 59\/3600\r\nlat1 == lat2\r\nlat1 - lat2\r\n<\/pre><pre class=\"codeoutput\">lat1 =\r\n         31.7\r\nlat2 =\r\n         31.7\r\nans =\r\n  logical\r\n   0\r\nans =\r\n  -3.5527e-15\r\n<\/pre><p>So a better idea might be to think about it geometrically. Two areas in the plane share a border if their intersection is a line or curve. And the new polygons in MATLAB have an <tt>intersect<\/tt> functionality! However, unfortunately, from a programming perspective, the intersection of two polygons is another polygon. And a polygon with zero area (which you get when the intersection is a line) is empty, so there's no obvious way to distinguish between states that share a boundary or not.<\/p><pre class=\"codeinput\">intersect(ma,nh)\r\nintersect(ma,me)\r\n<\/pre><pre class=\"codeoutput\">ans = \r\n  polyshape with properties:\r\n\r\n      Vertices: [0&times;2 double]\r\n    NumRegions: 0\r\n      NumHoles: 0\r\nans = \r\n  polyshape with properties:\r\n\r\n      Vertices: [0&times;2 double]\r\n    NumRegions: 0\r\n      NumHoles: 0\r\n<\/pre><p>Furthermore, this approach would still have been vulnerable to the finite precision problem.<\/p><h4>Solution: the fat of the land<a name=\"f61b2c99-b43b-4e09-9674-26d62c20a87c\"><\/a><\/h4><p>The new MATLAB polygons have a <tt>polybuffer<\/tt> function that gives you a way to \"fatten up\" a polygon by pushing the border out a given amount. Here's Massachusetts with a 0.1-degree (approx. 7-mile) border:<\/p><pre class=\"codeinput\">mafat = polybuffer(ma,0.1);\r\nplot([ma mafat])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_06.png\" alt=\"\"> <p>This gives a fairly simple, but robust, way to determine whether states are neighbors: put a small buffer around both and then see if there's any overlap.<\/p><pre class=\"codeinput\">pfat = polybuffer(p,1e-4); <span class=\"comment\">% add ~35 feet around every state<\/span>\r\nme = pfat(17);\r\nma = pfat(19);\r\nnh = pfat(27);\r\nintersect(ma,nh)\r\nintersect(ma,me)\r\n<\/pre><pre class=\"codeoutput\">ans = \r\n  polyshape with properties:\r\n\r\n      Vertices: [541&times;2 double]\r\n    NumRegions: 1\r\n      NumHoles: 0\r\nans = \r\n  polyshape with properties:\r\n\r\n      Vertices: [0&times;2 double]\r\n    NumRegions: 0\r\n      NumHoles: 0\r\n<\/pre><p>As you'd expect for a data type that represents polygons, there's an <tt>area<\/tt> function that does exactly what its name suggests:<\/p><pre class=\"codeinput\">area(intersect(ma,nh))\r\narea(intersect(ma,me))\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n   0.00035462\r\nans =\r\n     0\r\n<\/pre><p>So determining whether the padded states overlap should be simple: does the intersection have nonzero area? But the programmer in me was a bit worried: should I use a function, like <tt>area<\/tt>, that will come with some overhead? Would it be more efficient to directly query the properties of the polygon that results from the intersection?<\/p><p>However, there's another slight wrinkle that needs to be considered before worrying about implementation details. The <a href=\"https:\/\/en.wikipedia.org\/wiki\/Four_Corners\">\"Four Corners\"<\/a> area of the American Southwest is the one place where four states touch - Colorado, New Mexico, Arizona, and Utah. Does Colorado border Arizona? Does Utah border New Mexico? We can let philosophers and mathematicians argue over that. But unless an infinitely thin Loren is driving an infinitely thin car, the answer, for our purposes, is no. There's no way in reality to drive from Colorado to Arizona without going through Utah or New Mexico.<\/p><pre class=\"codeinput\">plot(pfat,<span class=\"string\">'LineStyle'<\/span>,<span class=\"string\">'none'<\/span>)\r\naxis([-109.055 -109.035 36.99 37.01])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_07.png\" alt=\"\"> <p>But if we look at the intersection between the \"fattened\" states, there will be a small overlap between Colorado and Arizona. Note that it will be <i>very<\/i> small - only slightly bigger than my apartment in graduate school, which should probably be considered a quantum amount of geographic area. Any real border between states will be at least an order of magnitude larger than this. So now we have a simple and robust way to determine whether or not states share a border:<\/p><div><ol><li>Add a small border to all states<\/li><li>Compute the area of the intersection of states<\/li><li>States share a border if that intersection is greater than a threshold<\/li><\/ol><\/div><pre class=\"codeinput\">area(intersect(ma,nh)) &gt; 1e-6\r\narea(intersect(ma,me)) &gt; 1e-6\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n  logical\r\n   1\r\nans =\r\n  logical\r\n   0\r\n<\/pre><p>Now I just need to apply that to every pair of states. I could probably try to get fancy with <tt>arrayfun<\/tt>, but this is a case where the simple <tt>for<\/tt>-loop approach might be a better idea. In fact, loops can help me avoid redundant calculations. If State A borders State B, then State B borders State A and I don't have to check that.<\/p><pre class=\"codeinput\">border = zeros(48);\r\n<span class=\"keyword\">for<\/span> k = 1:48\r\n    <span class=\"keyword\">for<\/span> j = (k+1):48 <span class=\"comment\">% Only need j for ks we haven't checked yet<\/span>\r\n        border(j,k) = area(intersect(pfat(j),pfat(k))) &gt; 1e-6;\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>The result of this is a lower-triangular matrix of connections. Now I can build a graph from the connectivity matrix and solve the problem. One of the neat graph functions is <tt>degree<\/tt>, which gives the number of connections to each node.<\/p><pre class=\"codeinput\">G = graph(border,<span class=\"string\">'lower'<\/span>);\r\nhistogram(degree(G))\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_08.png\" alt=\"\"> <p>So there's one degree-one node in the graph. That's the one we want. But which state does that correspond to? The names are stored in the <tt>Name<\/tt> field of the original structure array. I'll extract them into a cell array by using one of my favorite MATLAB tricks. When indexing produces multiple results (not just an array with multiple elements), those results come back in the form of a comma-separated list. That is, if <tt>s<\/tt> is a structure array (for example), <tt>s.Thing<\/tt> is equivalent to <tt>s(1).Thing, s(2).Thing, ..., s(end).Thing<\/tt>. What's great about having a comma-separated list is that that's the form MATLAB uses for concatenating elements into an array (i.e., <tt>[a,b,c]<\/tt> or <tt>{a,b,c}<\/tt>) and for passing inputs into a function (i.e., <tt>myfun(a,b,c)<\/tt>). Hence, I can extract all the <tt>Name<\/tt> field elements of the structure array <tt>states<\/tt> and concatenate them into a cell array with one simple line of MATLAB:<\/p><pre class=\"codeinput\">stnames = {states.Name};\r\n<\/pre><p>And now, finally, I can find out which state borders only one other state:<\/p><pre class=\"codeinput\">stnames{degree(G) == 1}\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n    'Maine'\r\n<\/pre><h4>Visualizing the results<a name=\"8e6d4ebf-e568-4117-a1fb-7be89984d444\"><\/a><\/h4><p>Having the state names, I can go a little further and visualize both the map and the connectivity graph. To start, it's possible to name the nodes in a graph:<\/p><pre class=\"codeinput\">G = graph(border,stnames,<span class=\"string\">'lower'<\/span>);\r\n<\/pre><p>This is useful when you visualize a graph because the nodes are automatically labeled:<\/p><pre class=\"codeinput\">plot(G)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_09.png\" alt=\"\"> <p>But this graph doesn't show the states in a physical layout. It's also possible to specify the locations of the nodes when plotting a graph, but for that I need a point in each state to use as the node location. Right in the middle of the state would be good... Guess what? There's a polygon function for that: <tt>centroid<\/tt>.<\/p><pre class=\"codeinput\">[x,y] = centroid(p);\r\nplot(G,<span class=\"string\">'XData'<\/span>,x,<span class=\"string\">'YData'<\/span>,y)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_10.png\" alt=\"\"> <p>That's looking pretty good. It would be nice to add the map for reference.<\/p><pre class=\"codeinput\">plot(p)\r\nhold <span class=\"string\">on<\/span>\r\nplot(G,<span class=\"string\">'XData'<\/span>,x,<span class=\"string\">'YData'<\/span>,y)\r\nhold <span class=\"string\">off<\/span>\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_11.png\" alt=\"\"> <p>Nice! Something cool to notice about that code: two uses of <tt>plot<\/tt>, both of which were on \"exotic\", non-numeric data types (one polygon array and one graph). This is something that our management team work hard on - making sure functionality is consistent, so that, for example, <tt>plot<\/tt> works naturally on different kinds of data. It's a feature of MATLAB I greatly appreciate. It means that I, as a MATLAB user, don't have to expend brain energy on which plotting function, or even which variant of <tt>plot<\/tt>, I have to use. I can just use <tt>plot<\/tt> and get on with what I'm actually trying to do.<\/p><h4>Where to next?<a name=\"f06aac2c-7f31-4aea-b5b8-f4ce0e01d8cd\"><\/a><\/h4><p>I solved the brain teaser, and even went a bit further to make a nice visualization of the states and their connections. But there's one thing about my visualization that's still bugging me: the states are colored with seven unique colors in the order that the states appear in the polygon array (which happens to be alphabetical by name). This means that neighboring states are sometimes given the same color. For example, note that Nebraska and the Dakotas are elments 25, 32, and 39, which means they all end up the same color. I happen to recall that there's <a href=\"https:\/\/en.wikipedia.org\/wiki\/Four_color_theorem\">a mathematical theorem<\/a> that says that any map can be colored with only four unique colors, such that no neighboring regions have the same color.<\/p><p>To me, that sounds like a second interesting task for polygons and graphs,... but a story for another time. In the meantime, <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=2456#respond\">let us know<\/a> in the comments if you can think of somewhere in your work where the new <tt>polyshape<\/tt> type might be useful.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_3f24d93ae0f34c7f8c4b480fb527092a() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='3f24d93ae0f34c7f8c4b480fb527092a ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 3f24d93ae0f34c7f8c4b480fb527092a';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2017 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_3f24d93ae0f34c7f8c4b480fb527092a()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2017b<br><\/p><\/div><!--\r\n3f24d93ae0f34c7f8c4b480fb527092a ##### SOURCE BEGIN #####\r\n%% Loren's Excellent Adventure: Maps, Graphs, and Polygons\r\n%\r\n% <https:\/\/www.mathworks.com\/products\/new_products\/latest_features.html\r\n% R2017b> was released recently. Today's guest blogger, Matt Tearle,\r\n% discovered a fun application of one of the\r\n% <https:\/\/www.mathworks.com\/help\/matlab\/release-notes.html many new\r\n% features> - solving a map-based puzzle.\r\n\r\n%% A traveling MathWorker problem\r\n% <http:\/\/www.cartalk.com\/ \"Car Talk\"> recently featured a\r\n% <http:\/\/www.cartalk.com\/puzzler\/bobos-dont-look-back-tour brain teaser\r\n% submitted by a former MathWorker>. The short, slightly modified version\r\n% is: Loren has to drive from Delaware to every state in the contiguous US,\r\n% visiting each state exactly once; Loren's boss - let's call him Jack -\r\n% offers to meet her in the last state, which turns out to be his\r\n% birthplace; but how can Jack know what the last state will be, without\r\n% waiting to hear Loren's itinerary?\r\n%\r\n% I'm about to give away the answer, so if you want to figure it out\r\n% yourself, do so now. It's OK, I'll wait.\r\n%\r\n% Done? Great. Hopefully you realized that there's one and only one state\r\n% that borders only one other state. That state has to be the beginning or\r\n% end of any path. (You can't visit any state more than once so once you\r\n% enter from the single neighboring state, you're stuck.) The beginning is\r\n% Delaware (which borders several other states), so this unique\r\n% \"degree-one\" state must be the end of Loren's epic journey.\r\n\r\n%% Enter MATLAB\r\n% As usual for me, the really interesting aspect of this problem is: can I\r\n% solve it in MATLAB? This is a classic \"traveling salesman\" problem, so\r\n% the graph capabilities recently added to MATLAB should be of use. But the\r\n% first problem is how to build the graph that represents the connections\r\n% between the states. It turns out that another new feature will help.\r\n%\r\n% A new variable type for representing 2-D polygons was introduced in\r\n% R2017b. The |polyshape| function takes vectors representing _x_ and _y_\r\n% coordinates and returns a polygon with those points as its vertices. It\r\n% so happens that Mapping Toolbox comes with a data file containing the\r\n% latitudes and longitudes of the borders of the US states - just what I\r\n% need. With the |shaperead| function I can import this data in the form of\r\n% a structure array, with fields for the name of the state and the\r\n% latitude\/longitude of the border:\r\nstates = shaperead('usastatehi.shp');\r\n% Remove Alaska & Hawaii, and Washington DC\r\nstates([2,11,51]) = [];\r\n\r\n%%\r\n% For a single state, it's straightforward to extract the arrays in the |X|\r\n% (longitude) and |Y| (latitude) fields and use these as the inputs to\r\n% |polyshape|:\r\nx = states(1).X;\r\ny = states(1).Y;\r\np = polyshape(x,y)\r\n\r\n%%\r\n% Now that the state is represented as a polygon, I can use the functions\r\n% that can work with polygon, such as plotting:\r\nplot(p)\r\n\r\n%%\r\n% But to really get the benefit of working with polygons, I'd like to get\r\n% an array of them, where each element is a polygon representing one\r\n% state. I could extract the _x_ and _y_ components for each state in a\r\n% loop, but any time I see myself about to apply an operation to each\r\n% element of an array, I think \"time for |arrayfun|!\".\r\nstruct2poly = @(s) polyshape(s.X,s.Y); % define my struct -> polyshape operation\r\np = arrayfun(struct2poly,states) % apply to all elements of the struct STATES\r\n\r\n%%\r\n% Because this is MATLAB, functions like |plot| should work on arrays of\r\n% polygons, just as well as scalar polygons, right?\r\nplot(p)\r\n\r\n%%\r\n% Isn't that nice? Each polygon is plotted as its own shaded shape. If we\r\n% zoom in, we can see the solution to the puzzle.\r\naxis([-75 -66 42 48])\r\n\r\n%% Problem: who is my neighbor?\r\n% By eye we can see how states share borders, but how can we get MATLAB to\r\n% determine this in an automated way? Let's start with the simplest\r\n% problem: given two states (polygons), can we determine if they share a\r\n% border or not?\r\nme = p(17);\r\nma = p(19);\r\nnh = p(27);\r\nsubplot(1,2,1)\r\nplot([ma nh])\r\ntitle('Shared border')\r\nsubplot(1,2,2)\r\nplot([ma me])\r\ntitle('No shared border')\r\n\r\n%%\r\n% The polygons have a |Vertices| property that contains the border\r\n% locations. So a simplistic approach would be to see if the two polygons\r\n% have any vertices in common. However, this has a lot of potential\r\n% problems. Firstly, two shapes can share a boundary without sharing\r\n% vertices:\r\na = polyshape([0 1 1 0],[0 0 3 3]);\r\nb = polyshape([1 2 2 1],[1 1 2 2]);\r\nsubplot(1,1,1)\r\nplot([a b])\r\nintersect(a.Vertices,b.Vertices,'rows')\r\n\r\n%%\r\n% Secondly, even if the survey points all align, there's always the\r\n% potential for floating-point precision issues: how is 31\u00c2\u00b041'59\"\r\n% represented?\r\nlat1 = 31.699722222222222\r\nlat2 = 31 + 41\/60 + 59\/3600\r\nlat1 == lat2\r\nlat1 - lat2\r\n\r\n%%\r\n% So a better idea might be to think about it geometrically. Two areas in\r\n% the plane share a border if their intersection is a line or curve. And\r\n% the new polygons in MATLAB have an |intersect| functionality! However,\r\n% unfortunately, from a programming perspective, the intersection of two\r\n% polygons is another polygon. And a polygon with zero area (which you get\r\n% when the intersection is a line) is empty, so there's no obvious way to\r\n% distinguish between states that share a boundary or not.\r\nintersect(ma,nh)\r\nintersect(ma,me)\r\n\r\n%%\r\n% Furthermore, this approach would still have been vulnerable to the finite\r\n% precision problem.\r\n\r\n%% Solution: the fat of the land\r\n% The new MATLAB polygons have a |polybuffer| function that gives you a way\r\n% to \"fatten up\" a polygon by pushing the border out a given amount. Here's\r\n% Massachusetts with a 0.1-degree (approx. 7-mile) border:\r\nmafat = polybuffer(ma,0.1);\r\nplot([ma mafat])\r\n\r\n%%\r\n% This gives a fairly simple, but robust, way to determine whether states\r\n% are neighbors: put a small buffer around both and then see if there's any\r\n% overlap.\r\npfat = polybuffer(p,1e-4); % add ~35 feet around every state\r\nme = pfat(17);\r\nma = pfat(19);\r\nnh = pfat(27);\r\nintersect(ma,nh)\r\nintersect(ma,me)\r\n\r\n%%\r\n% As you'd expect for a data type that represents polygons, there's an\r\n% |area| function that does exactly what its name suggests:\r\narea(intersect(ma,nh))\r\narea(intersect(ma,me))\r\n\r\n%%\r\n% So determining whether the padded states overlap should be simple: does\r\n% the intersection have nonzero area? But the programmer in me was a bit\r\n% worried: should I use a function, like |area|, that will come with some\r\n% overhead? Would it be more efficient to directly query the properties of\r\n% the polygon that results from the intersection?\r\n%\r\n% However, there's another slight wrinkle that needs to be considered\r\n% before worrying about implementation details. The\r\n% <https:\/\/en.wikipedia.org\/wiki\/Four_Corners \"Four Corners\"> area of the\r\n% American Southwest is the one place where four states touch - Colorado,\r\n% New Mexico, Arizona, and Utah. Does Colorado border Arizona? Does Utah\r\n% border New Mexico? We can let philosophers and mathematicians argue over\r\n% that. But unless an infinitely thin Loren is driving an infinitely thin\r\n% car, the answer, for our purposes, is no. There's no way in reality to\r\n% drive from Colorado to Arizona without going through Utah or New Mexico.\r\n%%\r\nplot(pfat,'LineStyle','none')\r\naxis([-109.055 -109.035 36.99 37.01])\r\n%%\r\n% But if we look at the intersection between the \"fattened\" states, there\r\n% will be a small overlap between Colorado and Arizona. Note that it will\r\n% be _very_ small - only slightly bigger than my apartment in graduate\r\n% school, which should probably be considered a quantum amount of\r\n% geographic area. Any real border between states will be at least an\r\n% order of magnitude larger than this. So now we have a simple and robust\r\n% way to determine whether or not states share a border:\r\n%%\r\n% # Add a small border to all states\r\n% # Compute the area of the intersection of states\r\n% # States share a border if that intersection is greater than a threshold\r\n%%\r\narea(intersect(ma,nh)) > 1e-6\r\narea(intersect(ma,me)) > 1e-6\r\n\r\n%%\r\n% Now I just need to apply that to every pair of states. I could probably\r\n% try to get fancy with |arrayfun|, but this is a case where the simple\r\n% |for|-loop approach might be a better idea. In fact, loops can help me\r\n% avoid redundant calculations. If State A borders State B, then State B\r\n% borders State A and I don't have to check that.\r\nborder = zeros(48);\r\nfor k = 1:48\r\n    for j = (k+1):48 % Only need j for ks we haven't checked yet\r\n        border(j,k) = area(intersect(pfat(j),pfat(k))) > 1e-6;\r\n    end\r\nend\r\n\r\n%%\r\n% The result of this is a lower-triangular matrix of connections. Now I can\r\n% build a graph from the connectivity matrix and solve the problem. One of\r\n% the neat graph functions is |degree|, which gives the number of\r\n% connections to each node.\r\nG = graph(border,'lower');\r\nhistogram(degree(G))\r\n\r\n%%\r\n% So there's one degree-one node in the graph. That's the one we want.\r\n% But which state does that correspond to? The names are stored in the\r\n% |Name| field of the original structure array. I'll extract them into a\r\n% cell array by using one of my favorite MATLAB tricks. When indexing\r\n% produces multiple results (not just an array with multiple elements),\r\n% those results come back in the form of a comma-separated list. That is,\r\n% if |s| is a structure array (for example), |s.Thing| is equivalent to\r\n% |s(1).Thing, s(2).Thing, ..., s(end).Thing|. What's great about having a\r\n% comma-separated list is that that's the form MATLAB uses for\r\n% concatenating elements into an array (i.e., |[a,b,c]| or |{a,b,c}|) and\r\n% for passing inputs into a function (i.e., |myfun(a,b,c)|). Hence, I can\r\n% extract all the |Name| field elements of the structure array |states| and\r\n% concatenate them into a cell array with one simple line of MATLAB:\r\nstnames = {states.Name};\r\n\r\n%%\r\n% And now, finally, I can find out which state borders only one other\r\n% state:\r\nstnames{degree(G) == 1}\r\n\r\n%% Visualizing the results\r\n% Having the state names, I can go a little further and visualize both the\r\n% map and the connectivity graph. To start, it's possible to name the nodes\r\n% in a graph:\r\nG = graph(border,stnames,'lower');\r\n\r\n%%\r\n% This is useful when you visualize a graph because the nodes are\r\n% automatically labeled:\r\nplot(G)\r\n\r\n%%\r\n% But this graph doesn't show the states in a physical layout. It's also\r\n% possible to specify the locations of the nodes when plotting a graph, but\r\n% for that I need a point in each state to use as the node location. Right\r\n% in the middle of the state would be good... Guess what? There's a polygon\r\n% function for that: |centroid|.\r\n[x,y] = centroid(p);\r\nplot(G,'XData',x,'YData',y)\r\n\r\n%%\r\n% That's looking pretty good. It would be nice to add the map for\r\n% reference.\r\nplot(p)\r\nhold on\r\nplot(G,'XData',x,'YData',y)\r\nhold off\r\n\r\n%%\r\n% Nice! Something cool to notice about that code: two uses of |plot|, both\r\n% of which were on \"exotic\", non-numeric data types (one polygon array and\r\n% one graph). This is something that our management team work hard on -\r\n% making sure functionality is consistent, so that, for example, |plot|\r\n% works naturally on different kinds of data. It's a feature of MATLAB I\r\n% greatly appreciate. It means that I, as a MATLAB user, don't have to\r\n% expend brain energy on which plotting function, or even which variant of\r\n% |plot|, I have to use. I can just use |plot| and get on with what I'm\r\n% actually trying to do.\r\n\r\n%% Where to next?\r\n% I solved the brain teaser, and even went a bit further to make a nice\r\n% visualization of the states and their connections. But there's one thing\r\n% about my visualization that's still bugging me: the states are colored\r\n% with seven unique colors in the order that the states appear in the\r\n% polygon array (which happens to be alphabetical by name). This means that\r\n% neighboring states are sometimes given the same color. For example, note\r\n% that Nebraska and the Dakotas are elments 25, 32, and 39, which means\r\n% they all end up the same color. I happen to recall that there's\r\n% <https:\/\/en.wikipedia.org\/wiki\/Four_color_theorem a mathematical theorem>\r\n% that says that any map can be colored with only four unique colors, such\r\n% that no neighboring regions have the same color.\r\n%\r\n% To me, that sounds like a second interesting task for polygons and\r\n% graphs,... but a story for another time. In the meantime, \r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=2456#respond let us know> in\r\n% the comments if you can think of somewhere in your work where the new\r\n% |polyshape| type might be useful.\r\n\r\n##### SOURCE END ##### 3f24d93ae0f34c7f8c4b480fb527092a\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/polygonmaps1_11.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p><a href=\"https:\/\/www.mathworks.com\/products\/new_products\/latest_features.html\">R2017b<\/a> was released recently. Today's guest blogger, Matt Tearle, discovered a fun application of one of the <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/release-notes.html\">many new features<\/a> - solving a map-based puzzle.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2017\/10\/12\/lorens-excellent-adventure-maps-graphs-and-polygons\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[57,15,41,6],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2456"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=2456"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2456\/revisions"}],"predecessor-version":[{"id":2458,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2456\/revisions\/2458"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=2456"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=2456"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=2456"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}