.st0{fill:#FFFFFF;}

Mapping Toolbox in MATLAB 

 August 22, 2020

By  Roshini

What does a Mapping toolbox do?

Imagine that you are given the task of displaying the population of countries. Wouldn’t it be more attractive when this data is displayed in a map, rather than a table? This is precisely what the mapping toolbox helps us with. Before seeing about how to create maps and globes, let us see what a Mapping toolbox does and what its uses are. Mapping Toolbox provides algorithms, functions, and app for analysing geographic data and creating map displays in MATLAB. You can import vector and raster data from a wide range of file formats and web map servers. Geospatial data can be combined with base map layers from multiple sources in a single map display.

Primary functions of the toolbox

In this blog, we shall see a few essential functions of this toolbox. The main topic that we will be concentrating in this blog would be performing various 2D and 3D maps and globe displays. We are going to see the following implementations in MATLAB:

  1. Creating Maps Using MAPSHOW
  2. Display a Rotating Globe
  3. Creating an Interactive globe
  4. Visualizing the COVID-19 cases in a map

1) Creating Maps Using MAPSHOW

As a beginner and very first step, we shall see how we can create different varieties of road maps using the mapshow function in a simple code. As the first step in every code, we would use the clc, clear, and close commands to clear all the variables and close all the opened images.

clc
clear
close all

Map 1: Roads - A Geographic Data Structure

As a first map, we would display a geographic data structure array with lines representing roads. We have an inbuilt data in MATLAB, which is the road data of a place named Concord. The road data is stored in a shapefile.

What is a shapefile?

The shapefile format is a geospatial vector data format for geographic information system (GIS) software. The shapefile format can spatially describe vector features: points, lines, and polygons, representing, for example, water wells, rivers, and lakes. Each item usually has attributes that describe it, such as name or temperature.

roads = shaperead('concord_roads.shp');
figure
mapshow(roads);
xlabel('easting in meters')
ylabel('northing in meters')

CreatingMap1

The above map has the roads of Concord plotted in a simple map display.

Map 2: Roads with Custom LineStyle

It is displaying the same map as the previous roadmap, but with a customized linestyle. Here we are using the colon for this purpose.

figure
mapshow('concord_roads.shp','LineStyle',':');
xlabel('easting in meters')
ylabel('northing in meters')

Creating Map 2

This time the map is the same as the first one except that the roads are displayed using colons. This can be customized according to your wish.

Map 3: Roads with SymbolSpec

In the third map, we shall display the road shape, and render using a SymbolSpec. To learn about the concord_roads.shp dataset, read its associated concord_roads.txt metadata file, which describes the attributes.


type concord_roads.txt

% Query the attributes in this roads file.
roads = shaperead('concord_roads.shp')

% Find out how many roads fall in each CLASS.
histcounts([roads.CLASS],'BinLimits',[1 7],'BinMethod','integer')

% Find out how many roads fall in each ADMIN_TYPE.
histcounts([roads.ADMIN_TYPE],'BinLimits',[0 3],'BinMethod','integer')

creatingMap3
CreatingMap3
CreatingMap3
CreatingMap3

After reading the concords text file, we can see that the Concord roads are divided based on 2 criteria:

  1. Class type – This is to show what kind of road it is such as a highway or a major street or a minor street, etc. There are 7 types of classes ranging from 1 to 7.

  2. Admin type - This again shows what kind of a road it is by showing which place they connect, such as a state road or local road, etc. There are 4 types of roads ranging from 0 to 3

We are taking the histogram of both types to see how many roads are there in each. From the histogram output regarding the details of the Concord roads, we notice that there are no roads in this file that are CLASS 1 or 7, and no roads in either ADMIN_TYPE 1 or 2.

Now we need to create a SymbolSpec to:

  1. Color local roads (ADMIN_TYPE=0) black.

  2. Color state roads (ADMIN_TYPE=3) red.

  3. Hide very minor roads (CLASS=5,6).


% Override default property of the SymbolSpec
roadspec = makesymbolspec('Line',...
{'Default', 'Color','green'}, ...
{'ADMIN_TYPE',0, 'Color','black'}, ...
{'ADMIN_TYPE',3, 'Color','red'},...
{'CLASS',6, 'Visible','off'},...
{'CLASS',5, 'Visible','off'},...
{'CLASS',[1 4], 'LineWidth',2});
figure
mapshow('boston_roads.shp','SymbolSpec',roadspec);
xlabel('easting in meters')
ylabel('northing in meters')

This time we have the roads displayed with different colors according to the class that they belong to.

Map 4: GeoTIFF Image of Boston

As a map four, we shall display the Boston GeoTIFF image.

What is a GeoTIFF image?

The TIFF imagery file format can be used to store and transfer digital satellite imagery, scanned aerial photos, elevation models, and scanned maps.

GeoTIFF - The TIFF file that has geographic (or cartographic) data embedded as tags within the TIFF file. The geographic data of a GeoTIFF file can be used to position the image in the correct location and geometry on the screen of a geographic information display.


figure
mapshow boston.tif
axis image manual off

S = shaperead('boston_placenames.shp');
surveyFeetPerMeter = unitsratio('sf','meter');
for k = 1:numel(S)
    S(k).X = surveyFeetPerMeter * S(k).X;
    S(k).Y = surveyFeetPerMeter * S(k).Y;
end

% Display the placenames.
text([S.X], [S.Y], {S.NAME}, 'Color', [0 0 0], 'BackgroundColor',[0.9 0.9 0],'Clipping','on');


CreatingMap4

Finally, we have the map display of Boston with its satellite image along with the place names displayed.

2) Display a Rotating Globe

As a second application of the mapping toolbox, we shall see how to animate to produce a rotating globe.

clc
clear
close all

Now we shall set up a Globe display. The view is from above the North Pole.

figure
axesm('globe');
gridm('GLineStyle','-','Gcolor',[.7 .8 .9],'Grid','on')

DisplayRotatingGlobe-2

In the above figure, the globe appears similar to how it will look when we see it from above the north pole. It looks like a 2-D circle, rather than a 3-D sphere.

Next, we shall spin the globe by one revolution.

set(gca,'Box','off','Projection','perspective')
spin

For this, we have used the spin.m function in MATLAB.

% spin.m: Rotates a view around the equator one revolution
% in 5-degree steps. Negative step makes it rotate normally
% (west-to-east).
for i=90:-5:-270
    view(i,23.5); % Earth's axis tilts by 23.5 degrees
    drawnow
end

DisplayRotatingGlobe

We have the skeletal structure of the 3-D globe in the above image.Add the base color of copper color to give a feeling of land and spin at a slower rate.

base = zeros(180,360);
baseR = georefcells([-90 90],[0 360],size(base));
copperColor = [0.62 0.38 0.24];
hs = geoshow(base,baseR,'FaceColor',copperColor);
setm(gca,'Galtitude',0.025);
axis vis3d
spin

DisplayRotatingGlobe

Show the Earth in space. Blacken the figure background, turn off the three axes, and spin again.

clmo(hs)

load topo
topo = topo / (earthRadius('km')* 20);
hs = meshm(topo,topolegend,size(topo),topo);
demcmap(topo)

set(gcf,'color','black');
axis off;
spin

DisplayRotatingGlobe

The above figure shows the globe with its geographic colors of green and blue.

Applying lighting, which shifts as the planet rotates and makes it brighter.

camlight right
lighting Gouraud;
material ([.7, .9, .8])

DisplayRotatingGlobe

The above picture shows a brighter globe and also with more clear terrains.

3) Creating an Interactive globe

In the third application topic, we will see how we can construct a map of major world cities enhanced with coastlines and terrain. We shall also see how to interactively pick a location and get the name and location of the nearest city.

clc
clear
close all

Set up a Map Axes Object and Render a Global Elevation Grid. Construct the axes.

figure
axesm bries
text(.8, -1.8, 'Briesemeister projection')
framem('FLineWidth',1)

Load and display a 1-by-1-degree elevation grid.

load topo
geoshow(topo, topolegend, 'DisplayType', 'texturemap')
%Improve the Terrain Display. Get a colormap appropriate for elevation.
demcmap(topo)
%Make it brighter
brighten(.5)
%Add Simplified Coastlines. Load global coastline coordinates
load coastlines
%Generalize the coastlines to 0.25-degree tolerance
[rlat, rlon] = reducem(coastlat,coastlon, 0.25);
%Plot the coastlines in brown.
geoshow(rlat, rlon, 'Color', [.6 .5 .2], 'LineWidth', 1.5)

Plot City Locations with Red Point Markers

%Read a shapefile containing names of cities worldwide and their coordinates
%in latitude and longitude
cities = shaperead('worldcities', 'UseGeoCoords', true);

%Extract the point latitudes and longitudes with extractfield, and add them to the map
lats = extractfield(cities,'Lat');
lons = extractfield(cities,'Lon');
geoshow(lats, lons,...
    'DisplayType', 'point',...
    'Marker', 'o',...
    'MarkerEdgeColor', 'r',...
    'MarkerFaceColor', 'r',...
    'MarkerSize', 3)
text(-2.8,-1.8,'Major World Cities')

Select Cities Interactively

runCitySelectionLoop = true; % Set to true to run optional city selection loop
if(runCitySelectionLoop)
    h1 = text(-2.8, 1.9, 'Click on a dot for its city name. Press ENTER to stop');
    h2 = text(-2.8, 1.7, '');
    h3 = text(-2.8, 1.5, 'City Coordinates.');
    while true
        [selected_lat,selected_lon] = inputm(1);
        if isempty(selected_lat)
            break% User typed ENTER
        end
        d = distance(lats, lons, selected_lat, selected_lon);
        k = find(d == min(d(:)),1);
        city = cities(k);
        geoshow(city.Lat, city.Lon, ...
            'DisplayType','point',...
            'Marker','o', ...
            'MarkerEdgeColor','k', ...
            'MarkerFaceColor','y', ...
            'MarkerSize', 3)
        h2.String = city.Name;
        h3.String = num2str([city.Lat, city.Lon],'%10.2f');
    end
    disp('End of input.')
end



As the output of this code, we can see how the map will be projected in a Briesemeister projection. Also, here I have placed the cursor near Delhi, and we can see the output name mentioned as New Delhi along with its latitude and longitude coordinates displayed as 28.63 and 77.21.

InteractiveGlobe

4) Visualizing the COVID-19 cases in a map

The entire world is now in a hazardous situation as we see the number of COVID-19 cases rising in huge numbers daily. As a final application of this blog, I felt it would be relevant if we can plot and visualize the Covid-19 cases in a map using the mapping toolbox.

Geobubble mapping

Geobubble mapping format will be used for the visualizing purpose, where the cases will be shown in the form of a bubble in each country. The size and color of the bubble would indicate the severity of the disease. Highly affected countries would have a more giant bubble in the map compared to a country with a low number of cases. This is the geobubble mapping format.

clc
clear
close all

Getting data of the countries

From the spreadsheet link, we can extract the data of country name, the number of COVID-19 active cases, date of updating this data, and the number of deaths in that country. A function is created for this purpose named covid19data.


function [table_overall_a] = covid19data()
warning("off")
a_tableinfo=webread('https://spreadsheets.google.com/feeds/list/1lwnfa-GlNRykWBL5y7tWpLxDoCfs8BvzWxFjeOZ1YJk/1/public/values?alt=json');
% a_tableinfo
for i = 1:1:size(a_tableinfo.feed.entry,1)
    Country_a(i) = string(a_tableinfo.feed.entry{i}.gsx_country.x_t);
    Updated_time_a(i) = string(a_tableinfo.feed.entry{i}.updated.x_t);
    Confirmed_Case_a(i) = str2double(string(a_tableinfo.feed.entry{i}.gsx_confirmedcases.x_t));
    Death_a(i) = str2double(string(a_tableinfo.feed.entry{i}.gsx_reporteddeaths.x_t));
end
table_overall_a = table(Country_a',Updated_time_a',Confirmed_Case_a',Death_a');
table_overall_a.Properties.VariableNames = {'Country' 'Updated' 'Case' 'Death'};
end

The webread function is used to read the webpage that has the data of the spreadsheet. The data is separated into four different columns using the for loop using the variables Country_a, Updated_time_a, Confirmed_Case_a and Death_a. From this function, we have the data for plotting ready.

Visualization of Geo Bubble mapping

We would first get the latitude and longitude data of countries from  https://developers.google.com/public-data/docs/canonical/countries_csv

The data from the location mentioned above is saved as an excel sheet named Geolocation.xlsx. This excel sheet has the data of the country's names along with the number of active cases, death, and the date of information updated. Every country name from this excel sheet is compared with the table_overall_a, and its index is obtained using find() command.

%% Getting the data of countries
[table_overall_a] = covid19data;

%% Visualization of Geo Bubble Mapping
%Get the latitude and longitude of countries from
%https://developers.google.com/public-data/docs/canonical/countries_csv
table_location = readtable('Geolocation.xlsx');
for i=1:1:height(table_overall_a)
    index=find(strcmp(strtrim(table_overall_a.Country(i)),strtrim(table_location.name)));
    if ~isempty(index)
        latitude(i) = table_location.latitude(index);
        longitude(i) = table_location.longitude(index);
    else
        latitude(i) = 40;
        longitude(i) = 150;
    end
end

The index is now used to get the latitude and longitude locations from table_overall_a.

A new table named geotable is created along with the latitude, longitude, and number of active cases. The data from this table is plotted using the geobubble() command and also set the size variable property to the number of COVID-19 cases.

geotable = table(latitude',longitude',table_overall_a.Death);
geotable.Properties.VariableNames = {'latitude' 'longitude' 'Case'};
figure()
gb = geobubble(geotable,'latitude','longitude','SizeVariable','Case','Title','Geobubble Map For Covid-19 Deaths');

Finally, the color of the bubbles must vary according to the range of cases. So, here we mention the range for low, medium, high, and super high.

gb.SourceTable.Severity = discretize(table_overall_a.Death,[0 10 50 10000 200000],'categorical', {'Low', 'Medium', 'High', 'Super High'});
gb.ColorVariable = 'Severity';

VisualiseCovidCasesInMap

Conclusion

In this blog, we came across a few essential functions and implementations of the mapping toolbox. We also saw one application of it relevant to today’s situation. The mapping toolbox can also be extended to various other applications such as Web mapping, Terrain and elevation analysis, Geometric Geodesy and map projections, Data representation and transformations, etc.

Did you find some helpful content from our video or article and now looking for its code, model, or application? You can purchase the specific Title, if available, and instantly get the download link.

Thank you for reading this blog. Do share this blog if you found it helpful. If you have any queries, post them in the comments or contact us by emailing your questions to [email protected]. Follow us on LinkedIn Facebook, and Subscribe to our YouTube Channel. If you find any bug or error on this or any other page on our website, please inform us & we will correct it.

If you are looking for free help, you can post your comment below & wait for any community member to respond, which is not guaranteed. You can book Expert Help, a paid service, and get assistance in your requirement. If your timeline allows, we recommend you book the Research Assistance plan. If you want to get trained in MATLAB or Simulink, you may join one of our training modules. 

If you are ready for the paid service, share your requirement with necessary attachments & inform us about any Service preference along with the timeline. Once evaluated, we will revert to you with more details and the next suggested step.

Education is our future. MATLAB is our feature. Happy MATLABing!

About the author 

Roshini

An Electronics and Communication engineering student. Passionate about learning and sharing with others.

{"email":"Email address invalid","url":"Website address invalid","required":"Required field missing"}

Connect with MATLAB Helper ®

Follow: YouTube Channel, LinkedIn Company, Facebook Page, Instagram Page

Join Community of MATLAB Enthusiasts: Facebook Group, Telegram, LinkedIn Group

Use Website Chat or WhatsApp at +91-8104622179

>