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:
- NOAA Storm Prediction Center Severe Weather GIS (SVRGIS) Page: http://www.spc.noaa.gov/gis/svrgis/
- Natural Earth: http://www.naturalearthdata.com/downloads/
- Free GIS Datasets: https://freegisdata.rtwilson.com
- 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.