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]; NEstates = shaperead('usastatelo', 'UseGeoCoords', true, ... 'BoundingBox', [MapLonLimit' MapLatLimit']);
Set up map axes with a projection suitable to display the New England states.
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.1345e+03
Create an autumn colormap for the six New England states, and then use the flipud command to invert the matrix.
fall = flipud(autumn(numel(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
