Browsed by
Author: Brandon

Mapping the World’s Flight Routes

Mapping the World’s Flight Routes

Why?

I find global flight path maps mesmerizing. In a single image these maps present the vastness of our interconnected civilization on a global scale. These maps are by no means rare. However, they can be a fun programing exercise to generate.

Where to get flight data?

OpenFlights.org has assembled several databases relate to air travel. Their database on air routes and airports provided the necessary information to determine flight paths. The database contains nearly 60,000 routes from 3200 airports.

Flight paths were colored based on the departure continent. To correlate the departure airport, with a depature continent, an ISO-3166 database was utilized.

The aforementioned datasets were assembled as MATLAB Tables for use by the plotting script.

The shapefile for the world map was obtained from ArcGIS.

How to calculate the flight path.

The great-circle is the shortest distance between any two points traveling along the surface of a sphere. The great-circle is defined as a circle centered on the centroid of the sphere and orientated such that the perimeter of the circle passes through the two locations. The arc defined by the great-circle provides the shortest travel route often used for long-distance air travel.

Algebraic formulas for calculating the great circle distance as well as latitude and longitude values of intermediate points along the arc can be found here: https://www.movable-type.co.uk/scripts/latlong.html

A MATLAB function for calculating the great circle distance as well as intermediate points along the route is included at the bottom of this article.

The result

The map below shows the results of this exercise.

 

Creating AVI Videos in MATLAB

Creating AVI Videos in MATLAB

This video will demonstrate how to make a simple AVI video file of a MATLAB plot.

The following functions will be essential for this tutorial:

  • VideoWriter – This function create the video file. The input to the function is the video name.
  • getFrame – This function captures the current state of the figure which is to be written as a frame to the video file.
  • writeVideo This function write the capture frame to the video file.

I will also persistent variables. Persistent variables are local variables within a function that are retained for subsequent function calls. Nested functions will also be employed. Nested functions can be a bit confusing at first to a novice. However, they are an excellent way to keep code clean. Nested functions also have access to the variables within their parent function and thus it is not necessary to pass them through as inputs.

Source Code for Simple Video

function Made_Video_of_Equation

%% Setup the figure size and frame rate
xPixels = 1280;
FPS = 30; % frames per second
Ratio = 16/9;
yPixels = xPixels/Ratio;
AVI_Name = 'Simple_Video.avi';

T_Span = 5;

%% Create the figure
f = figure();
set(f, 'Units', 'pixels', 'Resize', 'off',...
    'MenuBar', 'none', 'ToolBar', 'none');

% Figure location on the screen, can be larger than screen resolution
f.Position = [50 50 xPixels yPixels];
hold on;
xlim([0 5]);
ylim([-1 1]);
grid('on');
xlabel('Time (s)');
ylabel('Amplitude');

%% Now Let's Plot
zeta = 0.1;
omega = 2*pi;
alpha = omega*zeta;
gamma = omega*sqrt(1-zeta^2);

t = linspace(0,5,T_Span*FPS);
y = exp(-alpha*t).*(alpha*sin(gamma*t)./gamma + cos(gamma*t));

for Frame = 1:length(t)
    
    % Update objects within figure and plot
    UpdateFigure(Frame)
    
    % Grab the frame and add to avi
    if Frame ~= length(t)
        CreateAvi(false)
    else
        CreateAvi(true)
    end
    pause(1/FPS);
end
    function UpdateFigure(Frame)
        % Persistent variables are stored between function calls and are
        % accessible in subsequent calls
        persistent h
        if Frame == 1
            % First frame, create object handle
            h = plot(t(1:Frame), y(1:Frame), 'LineWidth', 4);
        else
            set(h, 'XData', t(1:Frame), 'YData', y(1:Frame));
            set(gca, 'FontSize',18,...
                'FontName', 'Times New Roman');
            set(findall(gca,'type','text'),'FontSize',18,...
                'FontName', 'Times New Roman');
        end
        
    end

    function CreateAvi(FinalFrame)
        persistent CreateAVIFile
        frame = getframe(f);
        % Create the video
        
        if isempty(CreateAVIFile)
            % Create and open file during first call
            CreateAVIFile = VideoWriter(AVI_Name, 'Motion JPEG AVI');
            CreateAVIFile.Quality = 100;
            CreateAVIFile.FrameRate = FPS;
            open(CreateAVIFile);
            writeVideo(CreateAVIFile,frame)
        elseif FinalFrame
            % Close avi file on final frame
            writeVideo(CreateAVIFile,frame)
            close(CreateAVIFile)
        else
            writeVideo(CreateAVIFile,frame)
            % Already created, do nothing
        end
    end

end

The type of videos one can create with MATLAB can be rather complex. The video below is a portion which was used to plot and track orbital debris in Earth orbit (debris not included in this snippet). The source code is included below the video file. Credit should be given to Ryan Gray for the snippet of code used to map the texture map to the globe. Keep in mind, due to the large size of the image used for the texture map, this video may take a while to render. Render time can be decreased by reducing the texture map image size.

Source Code for Globe Video
Globe_Files

Weather Forecast

Weather Forecast

Getting the Current Weather Conditions

First, we will start off by using the Weather API to get the current weather conditions. This data will be obtained in a json format and read with the webread function.

%% Get the current weather conditions
% Assemble the URL
% http://api.wunderground.com/api/{key}/conditions/q/CA/San_Francisco.json

CallURL_current = ['http://api.wunderground.com/api/', key,...
 '/conditions/q/', State, '/' City, '.json'];

% Grab the data from the API
Current_Data = webread(CallURL_current, options);

The webread function will return the data in a structure format. The current weather conditions are stored as strings and numbers, so some formatting is required:

%% Display the current information
disp(Current_Data.current_observation.display_location.full);
disp('Current Conditions');
disp([num2str(Current_Data.current_observation.temp_f, 3) char(176) 'F, ',...
 Current_Data.current_observation.relative_humidity ' Humidity, ',...
 Current_Data.current_observation.weather]);
disp(['Wind: ' Current_Data.current_observation.wind_string, char(10)]);

The previous code will write the current weather conditions to the display window in the following style:

Houghton, MI
Current Conditions
21.5°F, 88% Humidity, Snow
Wind: From the South at 7.0 MPH Gusting to 11.0 MPH

Getting the Forecast Conditions

Next, we will acquire forecasted weather conditions using the same API. The call url for the forecast will be slightly different from previously. This API will provide the weather forecast for the current plus next three days:

%% Get the forecast
% http://api.wunderground.com/api/{key}/forecast/q/CA/San_Francisco.json
CallURL_forecast = ['http://api.wunderground.com/api/', key,...
'/forecast/q/', State, '/' City, '.json'];

% Grab the data from the API
Forecast_Data = webread(CallURL_forecast, options);

disp('Forecast');
% Display Forecast
for ii = 1:8
disp(Forecast_Data.forecast.txt_forecast.forecastday(ii).title);
disp([strrep(Forecast_Data.forecast.txt_forecast.forecastday(ii).fcttext,...
'. ', char(10)) char(10)])
end

The previous code will provide the following output style:

Forecast
Wednesday
Occasional snow showers
High 27F
Winds SSW at 5 to 10 mph
Chance of snow 70%.

Wednesday Night
Snow this evening will transition to snow showers overnight
Low 26F
Winds W at 20 to 30 mph
Chance of snow 70%
About one inch of snow expected.

...

Saturday Night
Mainly cloudy with snow showers around before midnight
Low 14F
Winds WNW at 10 to 20 mph
Chance of snow 40%.

Current Radar Map

Now it is time to acquire the animated radar map. The map can be obtained in either a *.png, *.gif, or *.swf format. For this tutorial we will request the *.gif format so that we do not need to concern ourselves with making the animation. Matlab can play gif videos using the implay function. However, I find it easier to simply open the image with a web browser using the web function.

%% Get the radar data
CallURL_radar = ['http://api.wunderground.com/api/', key,...
 '/animatedradar/q/', State, '/' City,...
 '.gif?newmaps=1&timelabel=1&timelabel.y=10&num=5&delay=50',...
 '&width=600&height=480'];
web(CallURL_radar,'-new','-notoolbar');

The obtained radar for Houghton, MI on December 28th 2016 is:

 

Plotting the Forecast Conditions

Finally, using the forecast weather results, I will plot the high and low temperature as well as the icon describing the weather for each day. The forecast data structure contains a field for each day containing an icon name for each day. The Weather Underground API documentation contains a detailed list of the icon sets available and how to use them: https://www.wunderground.com/weather/api/d/docs?d=resources/icon-sets

I will not display the code to generate this plot.

Download the Code

I’ve uploaded the code for this work to the MATLAB file exchange. You can access the upload here: https://www.mathworks.com/matlabcentral/fileexchange/60935-get-current-weather-conditions

Shapefiles in MATLAB

Shapefiles in MATLAB

A shapefile is a standard geospacial data format for storing location, shape, and attributes of geospacial features. The mapping toolbox in MATLAB has tools for working with these datasets. In this tutorial we will work with some of this functionality.

Shapefiles are typically composed of very large datasets. Shapefiles are actually composed of multiple required files that are shared in a compressed *.zip format. The three required files are:

  • *.shp – The main file that stores the geometry
  • *.shx – Stores index information for the geometry
  • *.dbf – Stores attributes for each shape

More information about shapefiles can be found here.

Examples of shapefiles can be found at the following sources:

  1. NOAA  Storm Prediction Center Severe Weather GIS (SVRGIS) Page: http://www.spc.noaa.gov/gis/svrgis/
  2. Natural Earth: http://www.naturalearthdata.com/downloads/
  3. Free GIS Datasets: https://freegisdata.rtwilson.com
  4. US Energy Information Agency: https://www.eia.gov/maps/layer_info-m.php

For this tutorial, I will use the tornado (torn.zip) and states (states.zip) datasets from NOAA SVRGIS.

The following functions will be essential to this tutorial:

  • shaperead – Reads in the shapefiles
  • makesymbolspec – Defines structure for plotting attributes
  • mapshow – Enables the shape data to be displayed

First, we will start out by loading the shapefiles into MATLAB. We will use the shaperead function for this task. For large datasets, the shaperead function makes it easy to load only specific sections of the data based on from a certain value range. This functionality will not be necessary for this example.

% Load the shape files
Tornadoes = shaperead('torn.shp');
States = shaperead('states.shp');

For this tutorial, we will color the tornado paths based on the magnitude of the tornado. There are two approaches to do this. The first uses the matlab makesymbospec function. This function outputs a structure containing the symbolization specification:

ColorMap = colormap(jet(6));

TornSpec = makesymbolspec('Line',...
 {'mag',0, 'Color',ColorMap(1,:)},...
 {'mag',1, 'Color',ColorMap(2,:)},...
 {'mag',2, 'Color',ColorMap(3,:)},...
 {'mag',3, 'Color',ColorMap(4,:)},...
 {'mag',4, 'Color',ColorMap(5,:)},...
 {'mag',5, 'Color',ColorMap(6,:)});

For datasets containing a plethora of specifications, it may be desirable for assemble the data within a loop. This can be done for this example as follows:

% Set colormap for Fujita scale
ColorMap = colormap(jet(6));
TornSpec.ShapeType = 'Line';
for ii = 0:5
 TornSpec.Color{ii+1,1} = 'mag';
 TornSpec.Color{ii+1,2} = ii;
 TornSpec.Color{ii+1,3} = ColorMap(ii+1,:);
end

Now we will plot the the tornado paths overlaid with the boundaries of the states.

%% Plot the entire United States
% Top Plot - entire country
ax(1) = subplot(3,2,[1:4]);
h_states = mapshow(States, 'FaceColor', [1 1 1]);
h_torn = mapshow(ax(1), Tornadoes, 'SymbolSpec',TornSpec);
axis(ax(1),'equal');
ax(1).XLim = [-2600000,2400000];
ax(1).YLim = [-1900000,1500000];
title('Tornadoes (1950-2015)');

% Legend Information, we will simulate a legend by create patch lines and
% hiding them so they do not appear on the plot
for ii = 0:5
 v_patch = [0 0; 1 0; 1 1; 0 1];
 f_patch = [1 2 3 4];
 h_patch(ii+1) = patch('Faces',f_patch,'Vertices',v_patch,...
 'FaceColor',ColorMap(ii+1,:), 'DisplayName', ['F' num2str(ii)]);
 h_patch(ii+1).Visible = 'off';
end
h_legend = legend('toggle');
h_legend.Location = 'south';
h_legend.Orientation = 'horizontal';
axis off

We will also plot the number of tornadoes per year and the relative frequency of each tornado strength. I will not present the entire code to complete this, but I will give you enough to get started:

%% Add secondary plots
% Tornadoes Per Year
ax(2) = subplot(3,2,5);
Years = [Tornadoes.yr];

The results of this plot can be viewed below.


You may notice a significant increase in the number of tornadoes per year starting around 1990. According to the NOAA storm prediction center, this increase is a result of an increased ability to detect tornadoes resulting from the adoption of Doppler Radar network.

Using ODE45 to Solve a Differential Equation

Using ODE45 to Solve a Differential Equation

Introduction

For this tutorial, I will demonstrate how to use the ordinary differential equation solvers within MATLAB to numerically solve the equations of motion for a satellite orbiting Earth.

For two-body orbital mechanics, the equation of motion for an orbiting object relative to a much heavier central body is modeled as:

\ddot{\vec{r}}=-\frac{\mu}{r^{3}}\vec{r}

Where μ is the gravitational parameter of earth (398600 km3/s2)

The state-space representation of this equation becomes:

\left[ \begin{matrix} {\dot{\vec{r}}} \\ {\dot{\vec{v}}} \\ \end{matrix} \right]=\left[ \begin{matrix} {{0}_{3\times 3}} & {{\text{I}}_{3\times 3}} \\ {\mu }/{{{r}^{3}}}\; & {{0}_{3\times 3}} \\ \end{matrix} \right]\left[ \begin{matrix} {\vec{r}} \\ {\vec{v}} \\ \end{matrix} \right]

MATLAB has many ODE solvers available for the coder. The selection of the appropriate solver is dependent on the type of ODE you are solving and the desired accuracy. For this problem, we will use the ode45 solver which uses a Runge-Kutta iterative method to achieve 4th and 5th order accuracy. I recommend that students write their own Runge-Kutta function to better understand this algorithm prior to adopting that MATLAB internal function. The ode45 function within MATLAB uses the Dormand-Prince formulation.

To understand the input parameters for the ode45 function, type “doc ode45” and “doc odeset” in the MATLAB command window.

Now Let’s Get Started

For this problem, the equation of motion for the satellite will be coded as an anonymous function. Anonymous functions are convenient since they can be written within another script of function, but they can only contain one executable statement. For this problem, anonymous functions are an excellent choice to code the state-space equation. NOTE: Any variable referenced by an anonymous function that is not an input to the function, must be define before the function within the code. For example: mu.

% Define Constants
mu = 3.986012E5; % km^3/s^2

% Define the equation of motion
EoM = @(t, x) [zeros(3,3), eye(3); -(mu/norm(x(1:3,1))^3)*eye(3), zeros(3,3)]*x;

% Determine Orbital Period, this is just to determine how long to simulate for
E = norm(V0)^2/2 - mu/norm(R0);
a = -mu/(2*E);
n = sqrt(mu/a^3);
T = 2*pi/n;

% Now lets plot the results
options = odeset('MaxStep', 5); % Set maximum step size as 5 seconds
[time, State] = ode45(EoM, [0 T], X0, options);

AnimatedOrbit

 

How Did I Know How Long to Integrate For?

Simply knowing R and V at any point in time and the gravitational parameter of the central body, the period of the orbit can be determined.

First the total energy of the orbit needs to be calculated. The total energy is composed of the kinetic and potential energies and remains constant throughout the orbit:

E=\frac{\left | \vec{v} \right |^{2}}{2} - \frac{\mu}{\left | \vec{r} \right |}

The semi-major axis of the orbit can then be determined:

a=\frac{\mu}{2E}

Finally, from the semi-major axis, the orbital period becomes:

T=2\pi \sqrt{\frac{a^{3}}{\mu}}

NOTE: The orbital period can only be determined for circular and elliptical orbits (E < 0).

Removing Items for Plot Legends

Removing Items for Plot Legends

In this tutorial, an object will be removed from a plot legend with the use of the hasbehavior command. Often when creating plots, annotations or trivial lines may be included which can be distracting when included on a legend. For this tutorial, an RLC circuit is simulated and a shaded region is included ± 0.1 A to communicate a convergence region. We will remove this shaded region from the plot legend while keeping the object on the plot.

For the lines that we desire remain on the legend, we will assign the line name through he DisplayName property for the line handle. A DisplayName can also be assigned during plotting using the DisplayName property as follows: plot(x,y, ‘DisplayName’, ‘Desired Name’).

% Now lets assign names to each line
for ii = 1:length(zeta)
    h1(ii).DisplayName = ['\zeta = ' num2str(zeta(ii))];
end

Full_Legend

The first legend entry will be removed using the hasbehavior command. Note: the Shaded_Handle variable represents the line handle associated with handle assigned to the shaded object, which was assigned on line 26 of the full source code.

% Remove the first legend entry associated with the A2 handle
hasbehavior(Shaded_Handle,'legend',false);
legend('show');

Now, when we toggle the legend, the legend entry for the shaded region is removed and only the lines remain.

Truncated_Legend

% RLC Coefficients & Solution
C = 1;
L = 1;
R = .5;

t = linspace(0,20,100);
zeta = 0.4:.2:0.8;
I = zeros(length(zeta), length(t));
ii = 1;
for jj = 1:length(zeta)
    w0 = 1/sqrt(L*C);
    wd = w0*sqrt(1-zeta(jj)^2);
    B2 = 5;
    B1 = 0;
    I(ii,:) = B1*exp(-alpha*t).*cos(wd*t) +...
        B2*exp(-alpha*t).*sin(wd*t);
    ii = ii + 1;
end

% Create Figure
figure();
hold('on');
grid('on');
xlabel('Time (s)'); ylabel('Current (A)');
h1 = plot(t, I, 'LineWidth', 2);
Shaded_Handle = patch([0, 20, 20, 0], [0.1, 0.1, -0.1, -0.1],...
    'blue', 'LineStyle', 'none', 'FaceAlpha', 0.5);
uistack(Shaded_Handle, 'bottom'); % Move Shaded_Handle to bottom

% Now lets assign names to each line
for ii = 1:length(zeta)
    h1(ii).DisplayName = ['\zeta = ' num2str(zeta(ii))];
end

% Remove the first legend entry associated with the A2 handle
hasbehavior(Shaded_Handle,'legend',false);
legend('show');
Creating Animated Plots in MATLAB

Creating Animated Plots in MATLAB

This tutorial will demonstrate how to create animated plots using MATLAB. This will be demonstrated through the use of a Fourier approximation of a square wave. The infinite series representing the Fourier approximation of a square wave is:
 
x_{square}(t)= \frac{4}{\pi} \sum_{k=1}^\infty \frac{\sin\left(2\pi (2k-1) ft \right)}{(2k-1)}  
We will now create an animated GIF showing the first 20 terms in this Fourier approximation. For this plot, the only information that changes for each series in title and the y-data for the Fourier approximation. This can be handled efficiently by only updating that information for each iteration of the loop.

%% Define waveform properties
f = 1;          % Frequency (Hz)
t = 0:.001:2;   % Eval time (s)
Y_Fourier = zeros(size(t)); % Preallocate

%% Generate the frames
for k = 1:20;
    Y_Fourier = (4/pi)*sin(2*pi*(2*k-1)*f*t)/(2*k-1) + Y_Fourier;
    if k == 1
        % Only create the plot for 1st iteration, update wave for
        % subsequent
        figure();
        hold('on'); grid('on');
        ylim([-1.5 1.5]);
        set(gca, 'xtick', 0:0.5:2);
        h1 = plot(t, square(t*2*pi), 'LineWidth', 2);
        h2 = plot(t,Y_Fourier, 'LineWidth', 2);
        legend('Square Wave (1 Hz)', 'Fourier Approximation');
        xlabel('Time (s)');
        ylabel('Amplitude');
    else
        % Update y data
        set(h2, 'ydata', Y_Fourier)
    end
    title(['k = ' num2str(k)]);

Creating the frames

To create the frames for the animated GIF, we will save each plot of the Fourier approximation to a *.png file. It is possible to create an animation without first saving the frames to files by using the getframe() function within MATLAB. However, this method allows for easier control over the animation resolution as well as works in cases where the frames were created using other means, e.g. images obtained from a camera.

Each frame will be generated using the print() function. The necessary inputs for this process is the frame name, the image format, and the image resolution. For this example, the animation will be created with a resolution of 150 PPI. Keep in mind that this animation will be displayed on a screen and images with resolutions greater than the screen resolution will be displayed at the screen resolution.

    % Save as png with a resolution of 150 pixels per inch
    print(['Frame ' num2str(k)], '-dpng', '-r150');
end

Assembling the frames into an animated GIF image

Finally, now that the frames are created, we can assemble the into an animated GIF using the imwrite() command. Two of the necessary inputs for the imwrite command is the image data X and a color scale associated with that data included in map. The GIF image format only supports 8 bit images. If the image being read with imread() is greater than uint8, a conversion may be necessary.

GifName = 'SquareWaveAnimation.gif';
delay = 0.5;    % Delay between frames (s)
for ii = 1:20
    [A, ~] = imread(['Frame ' num2str(ii) '.png']);
    [X, map] = rgb2ind(A, 256);
    if ii == 1
        imwrite(X, map, GifName, 'gif', 'LoopCount', inf, 'DelayTime', delay)
    else
        imwrite(X, map, GifName, 'gif', 'WriteMode', 'append', 'DelayTime', delay)
    end
end

The result

The output from the above code becomes:

SquareWaveAnimation

Full code

%% Define waveform properties
f = 1;          % Frequency (Hz)
t = 0:.001:2;   % Eval time (s)
Y_Fourier = zeros(size(t)); % Preallocate

%% Generate the frames
for k = 1:20;
    Y_Fourier = (4/pi)*sin(2*pi*(2*k-1)*f*t)/(2*k-1) + Y_Fourier;
    if k == 1
        % Only create the plot for 1st iteration, update wave for
        % subsequent
        figure();
        hold('on'); grid('on');
        ylim([-1.5 1.5]);
        set(gca, 'xtick', 0:0.5:2);
        h1 = plot(t, square(t*2*pi), 'LineWidth', 2);
        h2 = plot(t,Y_Fourier, 'LineWidth', 2);
        legend('Square Wave (1 Hz)', 'Fourier Approximation');
        xlabel('Time (s)');
        ylabel('Amplitude');
    else
        % Update y data
        set(h2, 'ydata', Y_Fourier)
    end
    title(['k = ' num2str(k)]);
    % Save as png with a resolution of 150 pixels per inch
    print(['Frame ' num2str(k)], '-dpng', '-r150');
end

%% Now time to generate the animated gif
GifName = 'SquareWaveAnimation.gif';
delay = 0.5;    % Delay between frames (s)
for ii = 1:20
    [A, ~] = imread(['Frame ' num2str(ii) '.png']);
    [X, map] = rgb2ind(A, 256);
    if ii == 1
        imwrite(X, map, GifName, 'gif', 'LoopCount', inf, 'DelayTime', delay)
    else
        imwrite(X, map, GifName, 'gif', 'WriteMode', 'append', 'DelayTime', delay)
    end
end
Solving the Traveling Salesman Problem Using the Google Maps API

Solving the Traveling Salesman Problem Using the Google Maps API

The full source code for this problem will not be posted since my intent is not to write work that can easily be used in its entirety as a course project. If you are a course instructor and would like to use this as a demonstration, feel free to contact me with your request through your university email account.

Determining the distance between cities

Introduction

Using the Google Maps API, the driving distance, driving time, and GPS locations of each city can be determined. To do this, the Google Maps Distance Matrix and the Geocoding APIs was utilized. To request information from these APIs, one simply needs to generate a URL that contains a list of the information desired in a format specified by the API. When the URL is entered into a browser, the resulting webpage contains the desired information. A thorough description of how to use this API can be determined using the previous links.

An API key is not necessary to use the APIs; however, without a key, your IP address will be limited to 1000 requests per day (multiple users under the same IP will contribute to this limit). With an key, the limit is expanded to 2500 per key.

Now, consider the case where you have a 10 city traveling salesman problem. One could simply create a list of origin and destination cities which would result in a 10×10 matrix. The API could easily handle this request as long as the URL did not exceed the 2000 character limit. However, this is inefficient since the diagonals of this matrix are zeros and the matrix is symmetrical. I.e., the distance from A to B is the same as from B to A (obviously). Therefore, effective programming can be used here to minimize the number of API calls.

Reading the API data with MATLAB

The webread function within MATLAB makes it easy to grab and parse the data returned by the API.  This function can grab handle several return formats. For this example, we will use JSON format which the function will automatically recognize and the output variable will be a structure containing the desired data.

url=https://maps.googleapis.com/maps/api/distancematrix/json?origins=Madison+WI&amp;amp;destinations=Milwaukee+WI|Chicago+IL&amp;amp;units=imperial
DistanceData = webread(url);

To see the JSON file which MATLAB will parse, click here.

To better understand the format of the structure variable DistanceData run the previous code and study the format of the variable in the MATLAB workspace. Then try adding an additional origin city and see how the structure format changes.

For a six city problem, the following distance matrix and city center data was obtained using the API:

Driving Distance Between Cities (miles)

 MadisonDuluthClevelandChicagoHoustonDenver
Madison0329.38496.62147.41135.1962.42
Duluth329.380818.13468.911327.31065.8
Cleveland496.62818.130343.251296.81332.6
Chicago147.4468.91343.2501083.31003.7
Houston1135.11327.31296.81083.301033.6
Denver962.421065.81332.61003.71033.60

Travel Time Between Cities (Hours)

 MadisonDuluthClevelandChicagoHoustonDenver
Madison05.02947.60582.488917.28913.858
Duluth5.0294012.2367.118619.63415.061
Cleveland7.605812.23605.304219.44419.132
Chicago2.48897.11865.3042016.23414.271
Houston17.28919.63419.44416.234015.287
Denver13.85815.06119.13214.27115.2870
The background map was obtained through the use of the plot_google_map function written by Zohar Bar-Yehuda. Background image copyright Google Maps.
The background map was obtained through the use of the plot_google_map function written by Zohar Bar-Yehuda. Background image copyright Google Maps.