{"id":2205,"date":"2017-03-28T16:36:35","date_gmt":"2017-03-28T21:36:35","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=2205"},"modified":"2017-04-04T08:15:21","modified_gmt":"2017-04-04T13:15:21","slug":"find-that-signpost-using-optimization-and-the-google-maps-api-in-matlab-to-find-a-landmark","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2017\/03\/28\/find-that-signpost-using-optimization-and-the-google-maps-api-in-matlab-to-find-a-landmark\/","title":{"rendered":"Find That Signpost! Using optimization and the Google Maps API in MATLAB to find a landmark"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>Interesting problems are everywhere. Today's guest blogger, Matt Tearle, loves the chance to apply MATLAB to any intellectually stimulating puzzle. It's a fun way to learn and practice new features.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#32835939-5119-4a56-beb8-f3db32ef9c03\">A Postcard Becomes a Puzzle<\/a><\/li><li><a href=\"#26f8fe61-16ff-480e-b5df-bd8c6893ec87\">Let's Try It in MATLAB<\/a><\/li><li><a href=\"#c4215a12-d687-4448-86a3-419d90fffdc6\">Try an Optimization Approach?<\/a><\/li><li><a href=\"#f2a74a71-c8c1-4816-adf8-44f46ec1ab76\">What Does \"Distance\" Mean?<\/a><\/li><li><a href=\"#81762025-1514-42be-b00e-5b0653627c3d\">How Good is the Result?<\/a><\/li><li><a href=\"#61940804-83dc-4606-8b4e-6fde3642ecac\">What MATLAB Have You Discovered by \"Accident\"?<\/a><\/li><\/ul><\/div><h4>A Postcard Becomes a Puzzle<a name=\"32835939-5119-4a56-beb8-f3db32ef9c03\"><\/a><\/h4><p>My wife received this postcard in the mail:<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/postcard.png\" alt=\"\"> <\/p><p>This is a real signpost, directing drivers to various real towns in Maine.<\/p><p>Never one to pass up a teaching moment, I used this to explain to my daughter how you could (theoretically) figure out the location of the signpost if you knew where all these towns were, by finding the point that was the given distance from each town.<\/p><h4>Let's Try It in MATLAB<a name=\"26f8fe61-16ff-480e-b5df-bd8c6893ec87\"><\/a><\/h4><p>Why settle for abstract theory? Why not actually try it out? A little manual work with Google and Wikipedia and we (my long-suffering daughter and I) had the locations of the towns:<\/p><pre class=\"codeinput\">town = {<span class=\"string\">'Norway'<\/span>,<span class=\"string\">'Paris'<\/span>,<span class=\"string\">'Denmark'<\/span>,<span class=\"string\">'Naples'<\/span>,<span class=\"string\">'Sweden'<\/span>,<span class=\"string\">'Poland'<\/span>,<span class=\"string\">'Mexico'<\/span>,<span class=\"string\">'Peru'<\/span>,<span class=\"string\">'China'<\/span>};\r\nlon = -[70.540;70.501;70.804;70.610;70.823;70.393;70.543;70.406;69.517];\r\nlat =  [44.213;44.260;43.970;43.972;44.133;44.061;44.557;44.507;44.479];\r\n\r\nscatter(lon,lat)\r\ntext(lon,lat,town)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_01.png\" alt=\"\"> <p>Add a little context:<\/p><pre class=\"codeinput\">load <span class=\"string\">usastates.mat<\/span>\r\nmainelon = usastates(17).Lon;\r\nmainelat = usastates(17).Lat;\r\nplot(mainelon,mainelat)\r\nhold <span class=\"string\">on<\/span>\r\nscatter(lon,lat)\r\ntext(lon,lat,town)\r\nhold <span class=\"string\">off<\/span>\r\naxis <span class=\"string\">equal<\/span>\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_02.png\" alt=\"\"> <p>We can add the distances to each town.<\/p><pre class=\"codeinput\">r = [14;15;23;23;25;27;37;46;94];\r\n<\/pre><p>Now we need to find the location that is the given distance from each town. There is a circle of locations 14 miles from Norway. Another circle 15 miles from Paris, and so on. The signpost should be at the intersection of all those circles. So let's plot them and see.<\/p><p>But how do we get the circle of points? One way would be manual calculation, using the parametric equation for a circle. But a more accurate approach would take into account the fact that we live on a spherical planet. This means that degrees of latitude and longitude are not equally spaced at all locations. There are formulas for finding all the points equidistant from a point on the Earth, otherwise known as <i>small circles<\/i> (as opposed to <i>great circles<\/i>). But, as is so often the case with MATLAB, there's also <a href=\"https:\/\/www.mathworks.com\/help\/map\/ref\/scircle1.html\">a function to do it for you<\/a>, if you have the right toolbox. In this case, we're performing geographic calculations, so it's not surprising that Mapping Toolbox has what we need.<\/p><p>calculate lat\/lon of circles a given radius from a given center (using a distance measure of statute miles)<\/p><pre class=\"codeinput\">[clat,clon] = scircle1(lat,lon,r,[],earthRadius(<span class=\"string\">'sm'<\/span>));\r\n\r\n<span class=\"comment\">% plot results<\/span>\r\nplottowns(lat,lon,clat,clon,mainelat,mainelon)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_03.png\" alt=\"\"> <p>That plot doesn't look promising. There's no clear single intersection point.<\/p><pre class=\"codeinput\">xlim([-71 -70])\r\nylim([43.5 44.5])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_04.png\" alt=\"\"> <p>Note that <tt>plottowns<\/tt> is a function I wrote to display the town locations and the small circle around it (in the same color). In the past, I'd hesitate before cluttering up my hard drive with a function file for such a specific purpose. But R2016b delivered a long-requested feature: you can now add local functions to a script.<\/p><pre class=\"codeinput\">dbtype(<span class=\"string\">'mainesignpost.m'<\/span>,<span class=\"string\">'287:307'<\/span>)\r\n<\/pre><pre class=\"codeoutput\">\r\n287   \r\n288   function plottowns(lat,lon,clat,clon,maplat,maplon)\r\n289   % plot state border\r\n290   plot(maplon,maplat,'k')\r\n291   hold on\r\n292   \r\n293   % add town locations\r\n294   clr = parula(length(lon));  % make a set of colors for the lines\r\n295   scatter(lon,lat,[],clr)     % plot town locations, with given colors\r\n296   \r\n297   % add circles of distances\r\n298   l = plot(clon,clat);\r\n299   % change colors to match town markers\r\n300   for k = 1:length(l)\r\n301       l(k).Color = clr(k,:);\r\n302   end\r\n303   \r\n304   % adjust axes\r\n305   axis equal\r\n306   axis([-71.5 -68 43 46])\r\n307   hold off\r\n<\/pre><h4>Try an Optimization Approach?<a name=\"c4215a12-d687-4448-86a3-419d90fffdc6\"><\/a><\/h4><p>Maybe around (-70.7, 44.3) is a point that's close to correct for a few of the towns. And maybe it's not too far off for the rest? Given that the distances are rounded, and we can't know exactly where in the town they're measuring to, maybe it's expecting too much to look for a perfect solution. Instead, how about looking for the best fit? This becomes an optimization problem: find the location (lat\/lon) that minimizes the total error in all the distances.<\/p><p>So, given a location $x = (lat,lon)$, the objective function to minimize is<\/p><p>$$ d(x) = \\Sigma_k \\left( \\rho(x,y_k) - r_k \\right)^2$$<\/p><p>where $\\rho(x,y_k)$ is some distance function between the given point $x$ and the $k^{\\rm th}$ town location $y_k$.<\/p><p>Again, Mapping Toolbox <a href=\"https:\/\/www.mathworks.com\/help\/map\/ref\/distance.html\">provides the function for $\\rho$<\/a>:<\/p><pre class=\"codeinput\">distfun = @(x) distance(x(1),x(2),lat,lon);\r\n<\/pre><p>But if you read the documentation (and you always should!), you find that <tt>distance<\/tt> returns the distance in degrees of arc. Of course, Mapping Toolbox also provides <a href=\"https:\/\/www.mathworks.com\/help\/map\/ref\/deg2sm.html\">a function to convert from degrees to (statute) miles<\/a>:<\/p><pre class=\"codeinput\">distfun = @(x) deg2sm(distance(x(1),x(2),lat,lon),<span class=\"string\">'Earth'<\/span>);\r\n<\/pre><p>Now that I have $\\rho(x)$, it's easy to build the objective ($d(x)$) and run the optimization:<\/p><pre class=\"codeinput\">d = @(x) sum((distfun(x) - r).^2);\r\nx0 = [mean(lat) mean(lon)];  <span class=\"comment\">% Start in the geographic center of all the towns<\/span>\r\npostlocation = fmincon(d,x0)\r\n<\/pre><pre class=\"codeoutput\">\r\nLocal minimum possible. Constraints satisfied.\r\n\r\nfmincon stopped because the size of the current step is less than\r\nthe default value of the step size tolerance and constraints are \r\nsatisfied to within the default value of the constraint tolerance.\r\n\r\n\r\n\r\npostlocation =\r\n       44.164      -71.061\r\n<\/pre><p>Let's see the result.<\/p><pre class=\"codeinput\">plottowns(lat,lon,clat,clon,mainelat,mainelon)\r\n<span class=\"comment\">% add the predicted location<\/span>\r\nhold <span class=\"string\">on<\/span>\r\nplot(postlocation(2),postlocation(1),<span class=\"string\">'rp'<\/span>)\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\/mainesignpost_05.png\" alt=\"\"> <p>That's not right, surely. The postcard clearly said \"Maine signpost\". This location is in New Hampshire! Yes, there's some vaguery about where the town location is and rounding distances to the nearest mile, but that's not a big enough effect to see this level of inaccuracy.<\/p><p>I was heading down a rabbit-hole of optimization. Should I add constraints? Is it particularly sensitive to the starting location? Then I realized that I was missing something obvious and important.<\/p><h4>What Does \"Distance\" Mean?<a name=\"f2a74a71-c8c1-4816-adf8-44f46ec1ab76\"><\/a><\/h4><p>What does a distance on a signpost mean? Does it mean direct (\"as the crow flies\") distance? Or -- perhaps more likely, especially for a roadside sign -- driving distance?<\/p><p>That will change things considerably. (Small towns in Maine are not connected by straight highways.)<\/p><p>But how can we compute driving distance? Surely Mapping Toolbox can't do that. Well, no, but then I remembered seeing <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/60471-find-a-rendezvous-point\">an example<\/a> by another <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/profile\/authors\/1050816-will-campbell\">MathWorker, Will Campbell<\/a>, who used the Google Maps API to find driving time between two points. Adapting his code to use driving distance gives the local function <tt>drivedist<\/tt> that has the same interface as <tt>distance<\/tt>, but makes a web API call to Google Maps to determine the distance.<\/p><pre class=\"codeinput\">dbtype(<span class=\"string\">'mainesignpost.m'<\/span>,<span class=\"string\">'257:284'<\/span>)\r\n<\/pre><pre class=\"codeoutput\">\r\n257   %%\r\n258   function d = drivedist(lat1,lon1,lat2,lon2)\r\n259   % Google Maps API requires building a URL with approriate search terms\r\n260   base_url = 'https:\/\/maps.googleapis.com\/maps\/api\/directions\/json?';\r\n261   % Add starting location (lat1\/lon1) to URL\r\n262   origin = [num2str(lat1) ',' num2str(lon1)];\r\n263   base_url = [base_url 'origin=' origin '&amp;destination='];\r\n264   % Last part of the URL is your personal key. For security, I've saved this in a text file.\r\n265   keystr = ['&amp;key=' fileread('mykey.txt')];\r\n266   \r\n267   % Loop over destinations (lat2\/lon2) and get distances\r\n268   n = length(lat2);\r\n269   d = zeros(n,1);\r\n270   for k = 1:n\r\n271       % Add destination to URL. Finish URL with key\r\n272       destination = [num2str(lat2(k)) ',' num2str(lon2(k))];\r\n273       myurl = [base_url destination keystr];\r\n274       % Send request to Google Maps and unpack the json file returned\r\n275       dists = jsondecode(urlread(myurl));\r\n276       % Was the query successfully completed?\r\n277       if isequal(dists.status,'OK')\r\n278           % Yes. Good. Extract the distance value (which is in meters) and convert to miles\r\n279           d(k) = dists.routes(1).legs(1).distance.value\/1609.34;\r\n280       else\r\n281           % Nope. If it's a random glitch, let's just make that one value something big and hope for the best\r\n282           d(k) = 1000;\r\n283       end\r\n284   end\r\n<\/pre><p>A couple of notes about this function. Firstly, it demonstrates the benefits of sharing code. I would have had a hard time figuring out the use of the Google Maps API on my own. With Will's example as a starting point, I was able to adapt it to my purpose pretty easily. Secondly, new MATLAB features just keep making things easier (which means we can do more complex and interesting things!). Will used regular expressions to pull apart the JSON information returned by Google Maps. The <tt><a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/jsondecode.html\">jsondecode<\/a><\/tt> function was introduced in R2016b, which made it much easier to extract the driving distance I needed in my function.<\/p><p>Another MathWorker also pointed out that, with these newly aquired skills, I could have actually used the Google Maps API to automate finding the locations of the towns.<\/p><p>Now that I have a function for finding driving distance, it should be easy enough to simply switch out the distance function. But before handing this to <tt>fmincon<\/tt>, it's probably worth considering that the objective function may not vary smoothly with location. Location is continuous in latitude\/longitude, but driving distance depends on a discrete, finite set of roads. Rather than use a gradient-based solver, maybe a direct search method from Global Optimization Toolbox might work better.<\/p><pre class=\"codeinput\">distfun = @(x) drivedist(x(1),x(2),lat,lon);\r\nd = @(x) sum((distfun(x) - r).^2);\r\n\r\n[postlocation,dmin,~,optiminfo] = patternsearch(d,x0);\r\n<\/pre><p>Let's see the results<\/p><pre class=\"codeinput\">plottowns(lat,lon,clat,clon,mainelat,mainelon)\r\nhold <span class=\"string\">on<\/span>\r\nplot(postlocation(2),postlocation(1),<span class=\"string\">'rp'<\/span>)\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\/mainesignpost_06.png\" alt=\"\"> <p>That looks a lot more reasonable. Maybe I could try some different starting locations, but it's worth noting that this is expensive on function evaluations:<\/p><pre class=\"codeinput\">optiminfo.funccount\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n   195\r\n<\/pre><p>Each of those requires multiple calls to Google Maps, for a total of:<\/p><pre class=\"codeinput\">numel(lon)*optiminfo.funccount\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n        1755\r\n<\/pre><p>The Google Maps API allows only 2500 requests per day (unless you pay for it), so it's best not to go crazy.<\/p><p>Instead, let's see how good our solution is.<\/p><pre class=\"codeinput\">evaluateresults(distfun(postlocation),r,town)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_07.png\" alt=\"\"> <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_08.png\" alt=\"\"> <p>That's not bad! (And, yes, that's another local function I'm using to make the plots. I like this new feature!)<\/p><pre class=\"codeinput\">dbtype(<span class=\"string\">'mainesignpost.m'<\/span>,<span class=\"string\">'310:326'<\/span>)\r\n<\/pre><pre class=\"codeoutput\">\r\n310   \r\n311   function evaluateresults(ractual,rgiven,town)\r\n312   soln_vs_actual = [rgiven round(ractual)];\r\n313   figure\r\n314   bar(soln_vs_actual)\r\n315   xticklabels(town)\r\n316   xtickangle(60)\r\n317   ylabel('Distance (miles)')\r\n318   legend('Stated distance','Actual distance','Location','northwest')\r\n319   \r\n320   percerr = 100*abs(diff(soln_vs_actual,1,2)).\/rgiven;\r\n321   figure\r\n322   bar(percerr)\r\n323   ylim([0 100])\r\n324   xticklabels(town)\r\n325   xtickangle(60)\r\n326   ylabel('Percentage error')\r\n<\/pre><p>This function also uses some other new functions that I love: <tt>xticklabels<\/tt> and <tt>xtickangle<\/tt> for adding and rotating text labels on the axes.<\/p><h4>How Good is the Result?<a name=\"81762025-1514-42be-b00e-5b0653627c3d\"><\/a><\/h4><p>Of course, there's another way to find this signpost, also using Google: search for information on it! It turns out you can find it:<\/p><pre class=\"codeinput\">realpost = [44.244284,-70.785009];\r\nda = d(realpost);\r\nplottowns(lat,lon,clat,clon,mainelat,mainelon)\r\nhold <span class=\"string\">on<\/span>\r\nplot(postlocation(2),postlocation(1),<span class=\"string\">'rp'<\/span>)\r\nplot(realpost(2),realpost(1),<span class=\"string\">'mx'<\/span>)\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\/mainesignpost_09.png\" alt=\"\"> <p>That's looking good. Zoom in...<\/p><pre class=\"codeinput\">xlim([-71 -70])\r\nylim([43.5 44.5])\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_10.png\" alt=\"\"> <p>Wow! Nice job, MATLAB. The predicted location is just a couple of miles down the road from the real location:<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/compare.png\" alt=\"\"> <\/p><p>And, as a sanity check, how good could we actually do? What distances do we get from the distance function (with the town locations we're using) for the actual signpost location?<\/p><pre class=\"codeinput\">evaluateresults(distfun(realpost),r,town)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_11.png\" alt=\"\"> <img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_12.png\" alt=\"\"> <p>Interesting! It turns out that the signpost is wrong (at least according to Google Maps) -- Sweden, ME, is only about 13 miles away, not 25:<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/tosweden.png\" alt=\"\"> <\/p><p>Of course, it depends on your route, but the other option suggested by Google Maps is still only 18 miles. Maybe roads have changed since the signpost was made. But regardless, for the distances given and the current roads, the predicted location was about as good as it could possibly be.<\/p><pre class=\"codeinput\">disp(<span class=\"string\">'              Total error:'<\/span>)\r\ndisp(<span class=\"string\">'--------------------------------------'<\/span>)\r\ndisp(<span class=\"string\">'Predicted location  |  Actual location'<\/span>)\r\nfprintf(<span class=\"string\">'      %6.2f        |      %6.2f\\n'<\/span>,dmin,da)\r\n<\/pre><pre class=\"codeoutput\">              Total error:\r\n--------------------------------------\r\nPredicted location  |  Actual location\r\n      110.70        |      170.35\r\n<\/pre><h4>What MATLAB Have You Discovered by \"Accident\"?<a name=\"61940804-83dc-4606-8b4e-6fde3642ecac\"><\/a><\/h4><p>As you may have guessed, my poor daughter gave up on this quest long before I did. But I was learning new things. And new tricks in MATLAB are the best things to learn!<\/p><p>How about you? Have random problems or quixotic quests led you to new MATLAB knowledge? Share your favorite moments of MATLAB serendipity <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=2205#respond\">here<\/a>.<\/p><pre class=\"codeinput\"><span class=\"keyword\">function<\/span> d = drivedist(lat1,lon1,lat2,lon2)\r\n<span class=\"comment\">% Google Maps API requires building a URL with approriate search terms<\/span>\r\nbase_url = <span class=\"string\">'https:\/\/maps.googleapis.com\/maps\/api\/directions\/json?'<\/span>;\r\n<span class=\"comment\">% Add starting location (lat1\/lon1) to URL<\/span>\r\norigin = [num2str(lat1) <span class=\"string\">','<\/span> num2str(lon1)];\r\nbase_url = [base_url <span class=\"string\">'origin='<\/span> origin <span class=\"string\">'&amp;destination='<\/span>];\r\n<span class=\"comment\">% Last part of the URL is your personal key. For security, I've saved this in a text file.<\/span>\r\nkeystr = [<span class=\"string\">'&amp;key='<\/span> fileread(<span class=\"string\">'mykey.txt'<\/span>)];\r\n\r\n<span class=\"comment\">% Loop over destinations (lat2\/lon2) and get distances<\/span>\r\nn = length(lat2);\r\nd = zeros(n,1);\r\n<span class=\"keyword\">for<\/span> k = 1:n\r\n    <span class=\"comment\">% Add destination to URL. Finish URL with key<\/span>\r\n    destination = [num2str(lat2(k)) <span class=\"string\">','<\/span> num2str(lon2(k))];\r\n    myurl = [base_url destination keystr];\r\n    <span class=\"comment\">% Send request to Google Maps and unpack the json file returned<\/span>\r\n    dists = jsondecode(urlread(myurl));\r\n    <span class=\"comment\">% Was the query successfully completed?<\/span>\r\n    <span class=\"keyword\">if<\/span> isequal(dists.status,<span class=\"string\">'OK'<\/span>)\r\n        <span class=\"comment\">% Yes. Good. Extract the distance value (which is in meters) and convert to miles<\/span>\r\n        d(k) = dists.routes(1).legs(1).distance.value\/1609.34;\r\n    <span class=\"keyword\">else<\/span>\r\n        <span class=\"comment\">% Nope. If it's a random glitch, let's just make that one value something big and hope for the best<\/span>\r\n        d(k) = 1000;\r\n    <span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n\r\n\r\n<span class=\"keyword\">function<\/span> plottowns(lat,lon,clat,clon,maplat,maplon)\r\n<span class=\"comment\">% plot state border<\/span>\r\nplot(maplon,maplat,<span class=\"string\">'k'<\/span>)\r\nhold <span class=\"string\">on<\/span>\r\n\r\n<span class=\"comment\">% add town locations<\/span>\r\nclr = parula(length(lon));  <span class=\"comment\">% make a set of colors for the lines<\/span>\r\nscatter(lon,lat,[],clr)     <span class=\"comment\">% plot town locations, with given colors<\/span>\r\n\r\n<span class=\"comment\">% add circles of distances<\/span>\r\nl = plot(clon,clat);\r\n<span class=\"comment\">% change colors to match town markers<\/span>\r\n<span class=\"keyword\">for<\/span> k = 1:length(l)\r\n    l(k).Color = clr(k,:);\r\n<span class=\"keyword\">end<\/span>\r\n\r\n<span class=\"comment\">% adjust axes<\/span>\r\naxis <span class=\"string\">equal<\/span>\r\naxis([-71.5 -68 43 46])\r\nhold <span class=\"string\">off<\/span>\r\n<span class=\"keyword\">end<\/span>\r\n\r\n\r\n<span class=\"keyword\">function<\/span> evaluateresults(ractual,rgiven,town)\r\nsoln_vs_actual = [rgiven round(ractual)];\r\nfigure\r\nbar(soln_vs_actual)\r\nxticklabels(town)\r\nxtickangle(60)\r\nylabel(<span class=\"string\">'Distance (miles)'<\/span>)\r\nlegend(<span class=\"string\">'Stated distance'<\/span>,<span class=\"string\">'Actual distance'<\/span>,<span class=\"string\">'Location'<\/span>,<span class=\"string\">'northwest'<\/span>)\r\n\r\npercerr = 100*abs(diff(soln_vs_actual,1,2)).\/rgiven;\r\nfigure\r\nbar(percerr)\r\nylim([0 100])\r\nxticklabels(town)\r\nxtickangle(60)\r\nylabel(<span class=\"string\">'Percentage error'<\/span>)\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><pre class=\"codeoutput\">Optimization terminated: mesh size less than options.MeshTolerance.\r\n<\/pre><script language=\"JavaScript\"> <!-- \r\n    function grabCode_3540404c33f64bb69f0095e57434da29() {\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='3540404c33f64bb69f0095e57434da29 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 3540404c33f64bb69f0095e57434da29';\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_3540404c33f64bb69f0095e57434da29()\"><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; R2016b<br><\/p><\/div><!--\r\n3540404c33f64bb69f0095e57434da29 ##### SOURCE BEGIN #####\r\n%% Find That Signpost! Using optimization and the Google Maps API in MATLAB to find a landmark\r\n% Interesting problems are everywhere. Today's guest blogger, Matt Tearle, loves \r\n% the chance to apply MATLAB to any intellectually stimulating puzzle. It's a \r\n% fun way to learn and practice new features. \r\n%% A Postcard Becomes a Puzzle\r\n% My wife received this postcard in the mail:\r\n% \r\n% <<postcard.png>>\r\n% \r\n% This is a real signpost, directing drivers to various real towns in Maine.\r\n% \r\n% Never one to pass up a teaching moment, I used this to explain to my daughter \r\n% how you could (theoretically) figure out the location of the signpost if you \r\n% knew where all these towns were, by finding the point that was the given distance \r\n% from each town.\r\n%% Let's Try It in MATLAB\r\n% Why settle for abstract theory? Why not actually try it out? A little manual \r\n% work with Google and Wikipedia and we (my long-suffering daughter and I) had \r\n% the locations of the towns:\r\n%%\r\ntown = {'Norway','Paris','Denmark','Naples','Sweden','Poland','Mexico','Peru','China'};\r\nlon = -[70.540;70.501;70.804;70.610;70.823;70.393;70.543;70.406;69.517];\r\nlat =  [44.213;44.260;43.970;43.972;44.133;44.061;44.557;44.507;44.479];\r\n\r\nscatter(lon,lat)\r\ntext(lon,lat,town)\r\n%% \r\n% Add a little context:\r\n\r\nload usastates.mat\r\nmainelon = usastates(17).Lon;\r\nmainelat = usastates(17).Lat;\r\nplot(mainelon,mainelat)\r\nhold on\r\nscatter(lon,lat)\r\ntext(lon,lat,town)\r\nhold off\r\naxis equal\r\n%% \r\n% We can add the distances to each town.\r\n%%\r\nr = [14;15;23;23;25;27;37;46;94];\r\n%% \r\n% Now we need to find the location that is the given distance from each \r\n% town. There is a circle of locations 14 miles from Norway. Another circle 15 \r\n% miles from Paris, and so on. The signpost should be at the intersection of all \r\n% those circles. So let's plot them and see.\r\n% \r\n% But how do we get the circle of points? One way would be manual calculation, \r\n% using the parametric equation for a circle. But a more accurate approach would \r\n% take into account the fact that we live on a spherical planet. This means that \r\n% degrees of latitude and longitude are not equally spaced at all locations. There \r\n% are formulas for finding all the points equidistant from a point on the Earth, \r\n% otherwise known as _small circles_ (as opposed to _great circles_). But, as \r\n% is so often the case with MATLAB, there's also <https:\/\/www.mathworks.com\/help\/map\/ref\/scircle1.html \r\n% a function to do it for you>, if you have the right toolbox. In this case, we're \r\n% performing geographic calculations, so it's not surprising that Mapping Toolbox \r\n% has what we need.\r\n%%\r\n% calculate lat\/lon of circles a given radius from a given center\r\n% (using a distance measure of statute miles)\r\n[clat,clon] = scircle1(lat,lon,r,[],earthRadius('sm'));\r\n\r\n% plot results\r\nplottowns(lat,lon,clat,clon,mainelat,mainelon)\r\n%% \r\n% That plot doesn't look promising. There's no clear single intersection \r\n% point.\r\n%%\r\nxlim([-71 -70])\r\nylim([43.5 44.5])\r\n%% \r\n% Note that |plottowns| is a function I wrote to display the town locations \r\n% and the small circle around it (in the same color). In the past, I'd hesitate \r\n% before cluttering up my hard drive with a function file for such a specific \r\n% purpose. But R2016b delivered a long-requested feature: you can now add local \r\n% functions to a script.\r\n\r\ndbtype('mainesignpost.m','287:307')\r\n%% Try an Optimization Approach?\r\n% Maybe around (-70.7, 44.3) is a point that's close to correct for a few of \r\n% the towns. And maybe it's not too far off for the rest? Given that the distances \r\n% are rounded, and we can't know exactly where in the town they're measuring to, \r\n% maybe it's expecting too much to look for a perfect solution. Instead, how about \r\n% looking for the best fit? This becomes an optimization problem: find the location \r\n% (lat\/lon) that minimizes the total error in all the distances.\r\n% \r\n% So, given a location $x = (lat,lon)$, the objective function to minimize \r\n% is\r\n% \r\n% $$ d(x) = \\Sigma_k \\left( \\rho(x,y_k) - r_k \\right)^2$$\r\n% \r\n% where $\\rho(x,y_k)$ is some distance function between the given point $x$ \r\n% and the $k^{\\rm th}$ town location $y_k$.\r\n% \r\n% Again, Mapping Toolbox <https:\/\/www.mathworks.com\/help\/map\/ref\/distance.html \r\n% provides the function for $\\rho$>:\r\n%%\r\ndistfun = @(x) distance(x(1),x(2),lat,lon);\r\n%% \r\n% But if you read the documentation (and you always should!), you find that \r\n% |distance| returns the distance in degrees of arc. Of course, Mapping Toolbox \r\n% also provides <https:\/\/www.mathworks.com\/help\/map\/ref\/deg2sm.html a function \r\n% to convert from degrees to (statute) miles>:\r\n\r\ndistfun = @(x) deg2sm(distance(x(1),x(2),lat,lon),'Earth');\r\n%% \r\n% Now that I have $\\rho(x)$, it's easy to build the objective ($d(x)$) and \r\n% run the optimization:\r\n\r\nd = @(x) sum((distfun(x) - r).^2);\r\nx0 = [mean(lat) mean(lon)];  % Start in the geographic center of all the towns\r\npostlocation = fmincon(d,x0)\r\n%% \r\n% Let's see the result.\r\n\r\nplottowns(lat,lon,clat,clon,mainelat,mainelon)\r\n% add the predicted location\r\nhold on\r\nplot(postlocation(2),postlocation(1),'rp')\r\nhold off\r\n%% \r\n% That's not right, surely. The postcard clearly said \"Maine signpost\". \r\n% This location is in New Hampshire! Yes, there's some vaguery about where the \r\n% town location is and rounding distances to the nearest mile, but that's not \r\n% a big enough effect to see this level of inaccuracy.\r\n% \r\n% I was heading down a rabbit-hole of optimization. Should I add constraints? \r\n% Is it particularly sensitive to the starting location? Then I realized that \r\n% I was missing something obvious and important.\r\n%% What Does \"Distance\" Mean?\r\n% What does a distance on a signpost mean? Does it mean direct (\"as the crow \r\n% flies\") distance? Or REPLACE_WITH_DASH_DASH perhaps more likely, especially for a roadside sign \r\n% REPLACE_WITH_DASH_DASH driving distance?\r\n% \r\n% That will change things considerably. (Small towns in Maine are not connected \r\n% by straight highways.)\r\n% \r\n% But how can we compute driving distance? Surely Mapping Toolbox can't do \r\n% that. Well, no, but then I remembered seeing <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/60471-find-a-rendezvous-point \r\n% an example> by another <https:\/\/www.mathworks.com\/matlabcentral\/profile\/authors\/1050816-will-campbell \r\n% MathWorker, Will Campbell>, who used the Google Maps API to find driving time \r\n% between two points. Adapting his code to use driving distance gives the local \r\n% function |drivedist| that has the same interface as |distance|, but makes a \r\n% web API call to Google Maps to determine the distance.\r\n%%\r\ndbtype('mainesignpost.m','257:284')\r\n%% \r\n% A couple of notes about this function. Firstly, it demonstrates the benefits \r\n% of sharing code. I would have had a hard time figuring out the use of the Google \r\n% Maps API on my own. With Will's example as a starting point, I was able to adapt \r\n% it to my purpose pretty easily. Secondly, new MATLAB features just keep making \r\n% things easier (which means we can do more complex and interesting things!). \r\n% Will used regular expressions to pull apart the JSON information returned by \r\n% Google Maps. The |<https:\/\/www.mathworks.com\/help\/matlab\/ref\/jsondecode.html jsondecode>| \r\n% function was introduced in R2016b, which made it much easier to extract the \r\n% driving distance I needed in my function.\r\n% \r\n% Another MathWorker also pointed out that, with these newly aquired skills, \r\n% I could have actually used the Google Maps API to automate finding the locations \r\n% of the towns.\r\n% \r\n% Now that I have a function for finding driving distance, it should be easy \r\n% enough to simply switch out the distance function. But before handing this to \r\n% |fmincon|, it's probably worth considering that the objective function may not \r\n% vary smoothly with location. Location is continuous in latitude\/longitude, but \r\n% driving distance depends on a discrete, finite set of roads. Rather than use \r\n% a gradient-based solver, maybe a direct search method from Global Optimization \r\n% Toolbox might work better.\r\n%%\r\ndistfun = @(x) drivedist(x(1),x(2),lat,lon);\r\nd = @(x) sum((distfun(x) - r).^2);\r\n\r\n[postlocation,dmin,~,optiminfo] = patternsearch(d,x0);\r\n%% \r\n% Let's see the results\r\n\r\nplottowns(lat,lon,clat,clon,mainelat,mainelon)\r\nhold on\r\nplot(postlocation(2),postlocation(1),'rp')\r\nhold off\r\n%% \r\n% That looks a lot more reasonable. Maybe I could try some different starting \r\n% locations, but it's worth noting that this is expensive on function evaluations:\r\n\r\noptiminfo.funccount\r\n%% \r\n% Each of those requires multiple calls to Google Maps, for a total of:\r\n\r\nnumel(lon)*optiminfo.funccount\r\n%% \r\n% The Google Maps API allows only 2500 requests per day (unless you pay \r\n% for it), so it's best not to go crazy.\r\n% \r\n% Instead, let's see how good our solution is.\r\n%%\r\nevaluateresults(distfun(postlocation),r,town)\r\n%% \r\n% That's not bad! (And, yes, that's another local function I'm using to \r\n% make the plots. I like this new feature!)\r\n\r\ndbtype('mainesignpost.m','310:326')\r\n%% \r\n% This function also uses some other new functions that I love: |xticklabels| \r\n% and |xtickangle| for adding and rotating text labels on the axes.\r\n%% How Good is the Result?\r\n% Of course, there's another way to find this signpost, also using Google: search \r\n% for information on it! It turns out you can find it:\r\n%%\r\nrealpost = [44.244284,-70.785009];\r\nda = d(realpost);\r\nplottowns(lat,lon,clat,clon,mainelat,mainelon)\r\nhold on\r\nplot(postlocation(2),postlocation(1),'rp')\r\nplot(realpost(2),realpost(1),'mx')\r\nhold off\r\n%% \r\n% That's looking good. Zoom in...\r\n%%\r\nxlim([-71 -70])\r\nylim([43.5 44.5])\r\n%% \r\n% Wow! Nice job, MATLAB. The predicted location is just a couple of miles \r\n% down the road from the real location:\r\n% \r\n% <<compare.png>>\r\n% \r\n% And, as a sanity check, how good could we actually do? What distances do \r\n% we get from the distance function (with the town locations we're using) for \r\n% the actual signpost location?\r\n%%\r\nevaluateresults(distfun(realpost),r,town)\r\n%% \r\n% Interesting! It turns out that the signpost is wrong (at least according \r\n% to Google Maps) REPLACE_WITH_DASH_DASH Sweden, ME, is only about 13 miles away, not 25:\r\n% \r\n% <<tosweden.png>>\r\n% \r\n% Of course, it depends on your route, but the other option suggested by \r\n% Google Maps is still only 18 miles. Maybe roads have changed since the signpost \r\n% was made. But regardless, for the distances given and the current roads, the \r\n% predicted location was about as good as it could possibly be.\r\n%%\r\ndisp('              Total error:')\r\ndisp('REPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH')\r\ndisp('Predicted location  |  Actual location')\r\nfprintf('      %6.2f        |      %6.2f\\n',dmin,da)\r\n%% What MATLAB Have You Discovered by \"Accident\"?\r\n% As you may have guessed, my poor daughter gave up on this quest long before \r\n% I did. But I was learning new things. And new tricks in MATLAB are the best \r\n% things to learn!\r\n% \r\n% How about you? Have random problems or quixotic quests led you to new MATLAB \r\n% knowledge? Share your favorite moments of MATLAB serendipity\r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=2205#respond here>.\r\n\r\n%%\r\nfunction d = drivedist(lat1,lon1,lat2,lon2)\r\n% Google Maps API requires building a URL with approriate search terms\r\nbase_url = 'https:\/\/maps.googleapis.com\/maps\/api\/directions\/json?';\r\n% Add starting location (lat1\/lon1) to URL\r\norigin = [num2str(lat1) ',' num2str(lon1)];\r\nbase_url = [base_url 'origin=' origin '&destination='];\r\n% Last part of the URL is your personal key. For security, I've saved this in a text file.\r\nkeystr = ['&key=' fileread('mykey.txt')];\r\n\r\n% Loop over destinations (lat2\/lon2) and get distances\r\nn = length(lat2);\r\nd = zeros(n,1);\r\nfor k = 1:n\r\n    % Add destination to URL. Finish URL with key\r\n    destination = [num2str(lat2(k)) ',' num2str(lon2(k))];\r\n    myurl = [base_url destination keystr];\r\n    % Send request to Google Maps and unpack the json file returned\r\n    dists = jsondecode(urlread(myurl));\r\n    % Was the query successfully completed?\r\n    if isequal(dists.status,'OK')\r\n        % Yes. Good. Extract the distance value (which is in meters) and convert to miles\r\n        d(k) = dists.routes(1).legs(1).distance.value\/1609.34;\r\n    else\r\n        % Nope. If it's a random glitch, let's just make that one value something big and hope for the best\r\n        d(k) = 1000;\r\n    end\r\nend\r\nend\r\n\r\n\r\nfunction plottowns(lat,lon,clat,clon,maplat,maplon)\r\n% plot state border\r\nplot(maplon,maplat,'k')\r\nhold on\r\n\r\n% add town locations\r\nclr = parula(length(lon));  % make a set of colors for the lines\r\nscatter(lon,lat,[],clr)     % plot town locations, with given colors\r\n\r\n% add circles of distances\r\nl = plot(clon,clat);\r\n% change colors to match town markers\r\nfor k = 1:length(l)\r\n    l(k).Color = clr(k,:);\r\nend\r\n\r\n% adjust axes\r\naxis equal\r\naxis([-71.5 -68 43 46])\r\nhold off\r\nend\r\n\r\n\r\nfunction evaluateresults(ractual,rgiven,town)\r\nsoln_vs_actual = [rgiven round(ractual)];\r\nfigure\r\nbar(soln_vs_actual)\r\nxticklabels(town)\r\nxtickangle(60)\r\nylabel('Distance (miles)')\r\nlegend('Stated distance','Actual distance','Location','northwest')\r\n\r\npercerr = 100*abs(diff(soln_vs_actual,1,2)).\/rgiven;\r\nfigure\r\nbar(percerr)\r\nylim([0 100])\r\nxticklabels(town)\r\nxtickangle(60)\r\nylabel('Percentage error')\r\nend\r\n\r\n##### SOURCE END ##### 3540404c33f64bb69f0095e57434da29\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/2017\/mainesignpost_11.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>Interesting problems are everywhere. Today's guest blogger, Matt Tearle, loves the chance to apply MATLAB to any intellectually stimulating puzzle. It's a fun way to learn and practice new features.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2017\/03\/28\/find-that-signpost-using-optimization-and-the-google-maps-api-in-matlab-to-find-a-landmark\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[33,41,6],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2205"}],"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=2205"}],"version-history":[{"count":4,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2205\/revisions"}],"predecessor-version":[{"id":2277,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/2205\/revisions\/2277"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=2205"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=2205"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=2205"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}