Create Choropleth Map of Population Density

This example shows how to create a choropleth map of population density for the six New England states in the year 2000.

Import low-resolution U.S. state boundary polygons, setting the map limits for the New England region.

MapLatLimit = [41 48];
MapLonLimit = [-74 -66];

states = readgeotable('usastatelo.shp');
n = ["Maine" "New Hampshire" "Vermont" "Massachusetts" ...
    "Connecticut" "Rhode Island" "New York"];
rows = ismember(states.Name,n);
NEstates = states(rows,:);

Set up map axes with a projection suitable to display the New England states.

figure
axesm('MapProjection', 'eqaconic', 'MapParallels', [],...
  'MapLatLimit', MapLatLimit, 'MapLonLimit', MapLonLimit,...
  'GLineStyle', '-')
geoshow(NEstates, 'DisplayType', 'polygon', 'FaceColor','green')

Identify the maximum population density for New England states.

maxdensity = max([NEstates.PopDens2000])

maxdensity = 1.0032e+03

Create an autumn colormap for the six New England states, and then use the flipud command to invert the matrix.

fall = flipud(autumn(height(NEstates)));

Make a symbol specification structure, a symbolspec, that assigns an autumn color to each polygon according to the population density.

densityColors = makesymbolspec('Polygon', {'PopDens2000', ...
   [0 maxdensity], 'FaceColor', fall});

Display the map.

geoshow(NEstates, 'DisplayType', 'polygon', ...
   'SymbolSpec', densityColors)
title ({'Population Density in New England in 2000', ...
   'in Persons per Square Mile'})

Add a colorbar. You can also experiment with other colormaps.

caxis([0 maxdensity])
colormap(fall)
colorbar