# 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```