{"id":94,"date":"2021-04-30T12:19:45","date_gmt":"2021-04-30T17:19:45","guid":{"rendered":"https:\/\/wollen.org\/blog\/?p=94"},"modified":"2025-04-24T01:52:04","modified_gmt":"2025-04-24T06:52:04","slug":"working-with-the-new-2020-census-data","status":"publish","type":"post","link":"https:\/\/wollen.org\/blog\/2021\/04\/working-with-the-new-2020-census-data\/","title":{"rendered":"Working with the new 2020 Census data"},"content":{"rendered":"

Like a lot of data nerds I was excited to check out the new US Census apportionment data<\/a>. It only comes once per decade after all. I thought it would be a good chance to work on something I really enjoy: plotting maps.<\/p>\n

In this post we’ll take a simple statistic, state population growth from 2010 to 2020, and focus most of our energy on producing the plot. The full code will be available at the end of the post.<\/p>\n

Data sources:<\/strong><\/p>\n

    \n
  1. 2010 Census<\/a><\/li>\n
  2. 2020 Census<\/a><\/li>\n
  3. Shapefile<\/a><\/li>\n<\/ol>\n
    \n

    1. Prepare the data.<\/h4>\n

    Start with the imports. To produce this plot we’ll use GeoPandas<\/strong>, which is a package that (not surprisingly) looks and feels a lot like pandas<\/strong>. It keeps our data in a nice GeoDataFrame<\/em> object on which we can eventually call .plot()<\/code>.<\/p>\n

    import pandas as pd\r\nimport geopandas as gpd\r\nimport matplotlib.pyplot as plt\r\nfrom mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes\r\nfrom matplotlib import cm<\/pre>\n
    Since we’re plotting population change we’ll want the POPULATION2010<\/code> and POPULATION2020<\/code> columns in the same dataframe. Read both datasets and concatenate the dataframes together. Because axis=1<\/code> the dataframes are placed side-by-side rather than stacked one below the other. Fortunately our datasets have the same structure with state names in the same order, so we don’t have to worry about matching up rows.<\/p>\n
    df2010 = pd.read_csv(\"data\/2010.csv\")\r\ndf2020 = pd.read_csv(\"data\/2020.csv\")\r\n\r\ndf = pd.concat([df2010, df2020], axis=1)\r\ndf.set_index(\"STATE2010\", inplace=True)<\/pre>\n
    Create a column containing population change, then rate of change. This is the statistic we’ll use to color states on the map.<\/p>\n
    df.loc[:, \"CHANGE\"] = df[\"POPULATION2020\"] - df[\"POPULATION2010\"]\r\ndf.loc[:, \"CHANGE_RATE\"] = df[\"CHANGE\"] \/ df[\"POPULATION2010\"]<\/pre>\n
    At this point df.head(3)<\/code> looks like this:<\/p>\n
               POPULATION2010 STATE2020  POPULATION2020  CHANGE  CHANGE_RATE\r\nSTATE2010                                                               \r\nAlabama           4802982   Alabama         5030053  227071     0.047277\r\nAlaska             721523    Alaska          736081   14558     0.020177\r\nArizona           6412700   Arizona         7158923  746223     0.116366\r\n<\/pre>\n
    Since we’re plotting Alaska and Hawaii on separate inset axes we can’t use the built-in cmap<\/code> parameter. We’ll have to write our own function to generate colors. State colors will be placed in a COLOR<\/code> column that can be passed into plot()<\/code>. Higher population growth is represented by darker greens and if a population shrank from 2010 to 2020, the function returns a shade of red.<\/p>\n
    def get_color(rate, max_rate):\r\n    if rate >= 0:\r\n        greens_cmap = cm.get_cmap(\"Greens\", 100)\r\n        return greens_cmap(rate \/ max_rate)\r\n    else:\r\n        reds_cmap = cm.get_cmap(\"Reds\", 100)\r\n        return reds_cmap(abs(rate) \/ max_rate)\r\n\r\n\r\ndf.loc[:, \"COLOR\"] = df[\"CHANGE_RATE\"].apply(get_color, max_rate=df[\"CHANGE_RATE\"].max())\r\n<\/pre>\n
    Notice we haven’t yet created a GeoDataFrame<\/em>. Our shapefile uses the EPSG:4326<\/em> coordinate system, i.e. normal latitude\/longitude pairs. Although it isn’t necessary here, I’m in the habit of reprojecting my CRS to EPSG:3857<\/em>. Remember to filter out any shapefile elements that aren’t in df<\/code> (Puerto Rico and D.C.).<\/p>\n
    shapefile_path = \"shapefile\/cb_2018_us_state_20m.shp\"\r\ngdf = gpd.read_file(shapefile_path, epsg=4326)\r\ngdf = gdf.to_crs(epsg=3857)\r\ngdf = gdf[gdf[\"NAME\"].isin(df[\"STATE2020\"])<\/pre>\n
    Next we need to merge<\/code> the pandas<\/strong> dataframe, which contains color information, with the GeoDataFrame<\/em>. To accomplish this we’ll first rename a column so that GeoPandas<\/strong> knows which rows to match.<\/p>\n
    gdf.rename(columns={\"NAME\": \"STATE2020\"}, inplace=True)\r\ngdf = gdf.merge(df, on=\"STATE2020\")<\/pre>\n
    \n
    2. Plot the data on a map.<\/h4>\n
    Now for the interesting part: calling gdf.plot()<\/code>. On this axis we need to set x- and y-limits and limit our map to the continental United States. Alaska and Hawaii often present a problem when plotting a US map. The extreme west (east!) Aleutian Islands actually extend into the eastern hemisphere. So if we didn’t limit the window, GeoPandas<\/strong> would plot a zoomed-out figure containing all 360\u00b0 of longitude. But I don’t want to pick on them too much. We’ll add Alaska and Hawaii inset axes next.<\/p>\n
    fig, ax = plt.subplots(figsize=(12, 7.4))\r\nfig.subplots_adjust(left=0.0, right=1.0, bottom=0.0, top=0.955)\r\n\r\nborder_thickness = 0.4\r\nborder_color = \"#0F0F0F\"\r\n\r\nus_map = gdf.plot(ax=ax, color=gdf[\"COLOR\"], categorical=True,\r\n                  edgecolor=border_color, linewidth=border_thickness)\r\n\r\nus_map.set_xlim(-14.1e6, -7.24e6)\r\nus_map.set_ylim(2.5e6, 6.47e6)<\/pre>\n
    After a bit of tedious work the inset axes fall into place. This procedure is very similar to plotting the continental US above. We simply create two slices of gdf<\/code> restricted to Alaska and Hawaii respectively, then again call plot()<\/code>.<\/p>\n
    alaska_inset = zoomed_inset_axes(ax, zoom=0.21, bbox_to_anchor=(0.25, 0.29), bbox_transform=plt.gcf().transFigure)\r\nalaska_inset.set_xlim(-20.0e6, -14.36e6)\r\nalaska_inset.set_ylim(6.53e6, 11.75e6)\r\ngdf_alaska = gdf[gdf[\"STATE2020\"] == \"Alaska\"]\r\ngdf_alaska.plot(ax=alaska_inset, color=gdf[\"COLOR\"], categorical=True,\r\n                edgecolor=border_color, linewidth=border_thickness)\r\n\r\nhawaii_inset = zoomed_inset_axes(ax, zoom=1.05, bbox_to_anchor=(0.36, 0.22), bbox_transform=plt.gcf().transFigure)\r\nhawaii_inset.set_xlim(-17.87e6, -17.20e6)\r\nhawaii_inset.set_ylim(2.09e6, 2.57e6)\r\ngdf_hawaii = gdf[gdf[\"STATE2020\"] == \"Hawaii\"]\r\ngdf_hawaii.plot(ax=hawaii_inset, color=gdf[\"COLOR\"], categorical=True,\r\n                edgecolor=border_color, linewidth=border_thickness)<\/pre>\n
    Finally add a title and attribution text and turn off any axis labels that would clutter the map.<\/p>\n
    us_map.set_title(\"State Population Change  |  2010-2020\", {\"fontname\": \"Ubuntu Condensed\", \"fontsize\": 16})\r\nus_map.annotate(text=\"Data:  US Census Bureau\", xy=(-7.35e6, 2.54e6), size=11, family=\"Ubuntu Condensed\",\r\n                horizontalalignment=\"right\", verticalalignment=\"bottom\")\r\n\r\nus_map.set_axis_off()\r\nalaska_inset.set_axis_off()\r\nhawaii_inset.set_axis_off()\r\nplt.savefig(\"population_2010-2020.png\", facecolor=\"#FEFEFE\", dpi=300)<\/pre>\n
    I usually suggest saving maps in vectorized SVG<\/em> format rather than PNG<\/em> or JPG<\/em> whenever possible. This avoids blurriness when zooming in. If you do output a raster image consider increasing dpi<\/code> from the default of 100<\/p>\n
    The output:<\/strong><\/p>\n
    \"\"<\/a><\/p>\n
    \n
    TODO:<\/strong><\/p>\n