spmap: Display Alaska and Hawaii next to contiguous United States

[Jacqueline Fiore]

Attachments provided to compare maps that exclude and include Alaska and Hawaii.


I am having trouble using spmap to create a map that displays Alaska and Hawaii next to the contiguous United States. When I execute my do file that excludes Alaska and Hawaii, the map is not distorted. However, when I revise the code to include Alaska and Hawaii, the map that displays all states looks very distorted.

I provided my Stata code that was used to generate the map without Alaska and Hawaii. How can I revise this code to generate a map that is not distorted and includes Alaska and Hawaii?

************************************************************************;
*Read in the data;
*Note: the master dataset excludes Puerto Rico;
************************************************************************;

tempfile tf;

use `moddata'master,clear;

collapse
(max)
ME
controls_noprior
controls_prior
,by(state);

gen     coverage=1 if ME==1;
replace coverage=2 if controls_noprior==1;
replace coverage=3 if controls_prior==1;

rename state stname;

save `tf', replace;

************************************************************************;
*TRANSLATE SHAPE FILES INTO .DTA;
************************************************************************;

cd `mapdata';

shp2dta using gz_2010_us_040_00_20m, database(states-d) coordinates(states-c) genid(_ID) replace;

use states-c;
    joinby _ID using states-d;
    drop if NAME=="Hawaii" | NAME=="Alaska" | NAME=="Puerto Rico";
    rename NAME stname;
save states-c, replace;

use states-d, clear;

rename NAME stname;

joinby stname using `tf', unmatched(both);
    tab stname if _merge!=3;
    drop if _merge!=3;
    drop _merge;

drop if stname=="Hawaii" | stname=="Alaska";
    
    
************************************************************************;
*SPMAP;
************************************************************************;

local mapvar="coverage";
spmap `mapvar' using states-c,
    freestyle aspect(0.5, placement(north)) id(_ID) clmethod(unique) fcolor(gs5 gs16 gs12) 
    yscale(off) xscale(off) ylabel(minmax, nogrid nolabels) xlabel(minmax, nogrid nolabels)
    plotregion(color(white)) graphregion(color(white))  
    legend(size(*1.45) ring(1) position(6)
    label(2 "Medicaid Expansion States (N=30)")
    label(3 "Non-Expansion States with No Prior Coverage (N=15)")
    label(4 "Non-Expansion States with Prior Coverage (N=6)")             
    );
    graph export "`latex'map_v3.png", replace;

[Scott Merryman]

The Aleutian Islands cross the international date line and wrap around the map. One way to deal with is to shrink Alaska and move Alaska and Hawaii to the SW. Below is an example:

cd "C:\Users\scott.merryman\Desktop\mapexample"
copy http://www2.census.gov/geo/tiger/GENZ2015/shp/cb_2015_us_state_500k.zip  cb_2015_us_state_500k.zip
unzipfile cb_2015_us_state_500k.zip, replace 
shp2dta using  cb_2015_us_state_500k, data("us_data")  coor("us_coordinates") /// 
 genid(id) gencentroids(c) replace
use us_data, clear
quietly destring STATEFP, generate(fips)
drop if fips > 56
save usstates,replace

//Determine IDs for Hawaii and Alaska
use usstates,clear
l fips id  if fips == 2 | fips == 15

// Alaska is reduced in size by 60% and then both Alaska and Hawaii are re-positioned under the southwestern states.
//Drop northwest islands of Hawaii and western part of Aleutian Islands 
use us_coordinates,clear

gen order = _n
drop if _X < -165 & _X != . &  _ID == 50
replace _X = _X  + 55  if  _X != .  &  _ID == 50
replace _Y = _Y  + 4  if _Y != .  &  _ID == 50

replace _X = _X*.4  -55 if  _X !=  . &  _ID == 15
replace _Y = _Y*.4  + 1 if _Y != .  & _ID == 15
drop if _X > -10 & _X != . & _ID == 15
sort order 
sort _ID
drop order

save us_coordinates2,replace

use usstates,clear
spmap using us_coordinates2  , id(id)

[Jacqueline Fiore]

Thank you! I did not realize the Aleutian Islands were the cause of the distorted map. I will apply your suggested changes this evening and post an update to confirm the fix worked.

[Robert Picard]

Mapping raw geographic coordinates (lat/lon) without projection is like looking at standard definition TV in stretch mode. It's not flattering to the actors.

You can use geo2xy (from SSC) to project such coordinates. To install geo2xy and download its supporting datasets in the current directory, type in Stata's Command window:

ssc install geo2xy
net get geo2xy

By default, geo2xy uses a Web Mercator projection, the same used by most web mapping applications (Google Maps, Bing, etc.). In the following example, I use an Albers projection, a projection well suited for maps showing the conterminous United States. Both Alaska and Hawaii are projected separately and once in (x,y) coordinates, moving them is just a matter of scale and offset. The projection is in meters, which explains the large offsets. A bit of trial and error is required to find a good composition. I would probably add a frame around each to make clear that these are insets.

use "geo2xy_us_coor.dta", clear

* drop Puerto Rico and other territories
drop if inlist(_ID, 55,54,32,28,29)

* Alaska
drop if _X > 0 & !mi(_X) & _ID == 37  // land in the Eastern hemisphere
geo2xy _Y _X if _ID == 37, replace proj(albers)
replace _Y = _Y / 3 + 700000 if _ID == 37
replace _X = _X / 3 - 1700000 if _ID == 37

* Hawaii
drop if _X < -160 & !mi(_X) & _ID == 40  // small islands to the west
geo2xy _Y _X if _ID == 40, replace proj(albers)
replace _Y = _Y / 2 + 1500000 if _ID == 40
replace _X = _X / 2 - 800000 if _ID == 40

* contiguous US
geo2xy _Y _X if !inlist(_ID, 37,40), replace proj(albers)

save "xy_coor.dta", replace

use "geo2xy_us_data.dta",clear
spmap using "xy_coor.dta", id(_ID)

The resulting map looks like:

[Jacqueline Fiore]

Hi Scott,

Thank you for providing the code and explanation. I revised my code and the new map looks great.

Kind regards,
Jacqueline

[Jacqueline Fiore]

Hi Robert,
Thank you for the helpful feedback. I will update my code with the geo2xy command.
Many thanks!
Jacqueline

[Felix Quinones]

Hello All,

I know that this post was already answered, but I want to do exactly the same, but only for Puerto Rico and its 78 municipalities.
I'm being trying to make the map using spmap, but every information that I found exclude Puerto Rico. Is there a way to create this maps for Puerto Rico?

Thanks,
Felix

[Robert Picard]

The problem here is to find an appropriate shapefile that contains the boundaries for the 78 municipalities. I must admit that this has been harder than I thought it should be. I found one that seems to work here (click on the PRI_adm1.zip download button far to the right). The license is "Humanitarian use only". Once unzipped and with Stata's current directory set to the one that contains the files, the following code will create a map of the boundaries of Puerto Rico's 78 municipalities:

* Source: https://data.humdata.org/dataset/peurto-rico-admin0-and-admin1
* PRI_adm1.zip (132.2K)
* Puerto Rico administrative level 1 (municipio) boundaries.
clear
shp2dta using "PRI_adm1.shp", data("pr_data.dta")  coor("pr_coor.dta") replace

* apply a Google Maps projection
use "pr_coor.dta", clear
geo2xy _Y _X, replace
save "pr_coor_XY.dta", replace

* create the map
use "pr_data", clear
spmap using "pr_coor_XY.dta", id(_ID)
graph export "pr_coor_XY.png", width(1000) replace

and the resulting map:

[Thomas Wilk]

Hi Robert, I appreciate you offering the code, however I'm having trouble reproducing the Albers projection map as you've posted. My map looks quite different. I'm not sure what could be different, I literally copied your code into Stata. Any insight would be appreciated.

[Robert Picard]

The geo2xy package has been updated since I posted #4. If I recall correctly, I switched to a US shapefile with less detail and therefore the polygon identifiers (_ID) do not match the ones used in #4. There's a new and improved example in the "Fun with Maps" section of the help file. Once you download the ancillary datasets, you run:

use "geo2xy_us_coor.dta", clear

// flip longitudes to reconnect Hawaii and Alaska
replace _X = cond(_X > 0, _X - 180, _X + 180) if inlist(_ID, 14, 42)

// Alaska - USGS recommends standard parallels of 55 and 65 north
sum _X if _ID == 14
local midlon = (r(min) + r(max)) / 2
geo2xy _Y _X if _ID == 14, replace ///
        proj(albers, 6378137 298.257223563 55 65 0 `midlon')
replace _Y = _Y / 3 + 800000 if _ID == 14
replace _X = _X / 3 - 1700000 if _ID == 14

// Hawaii - USGS recommends standard parallels of 8 and 18 north
sum _X if _ID == 42
local midlon = (r(min) + r(max)) / 2
geo2xy _Y _X if _ID == 42, replace ///
        proj(albers, 6378137 298.257223563 8 18 0 `midlon')
replace _Y = _Y / 1.2 + 850000 if _ID == 42
replace _X = _X / 1.2 - 800000 if _ID == 42

// Puerto Rico
geo2xy _Y _X if _ID == 39, replace proj(albers)
replace _Y = _Y + 500000 if _ID == 39
replace _X = _X + 2000000 if _ID == 39

// 48 states - USGS recommends standard parallels of 29.5 and 45.5 north
sum _X if !inlist(_ID, 14, 42, 39)
local midlon = (r(min) + r(max)) / 2
geo2xy _Y _X if !inlist(_ID, 14, 42, 39), replace ///
        proj(albers, 6378137 298.257223563 29.5 45.5 0 `midlon')

save "xy_coor.dta", replace

use "geo2xy_us_data.dta",clear
spmap using "xy_coor.dta", id(_ID)  name(fun_composite, replace)

and the resulting map looks like:

[Thomas Wilk]

Thank you Robert! Is there a way to tell what states and territories match up with the ID? Also, do you have a similar map for county and zip code shapes? Thank you!

[Robert Picard]

You convert a shapefile for use with Stata using shp2dta (from SSC) or with the new spshape2dta command (Stata 15 or higher). These create two Stata datasets, one for coordinates that you use to draw the polygons and the other for the "database", one observation per _ID. Here's a sample of the matching database for the data used in #10:

. use "geo2xy_us_data.dta", clear

. list STATEFP STUSPS NAME _ID in 1/5

     +-------------------------------------+
     | STATEFP   STUSPS         NAME   _ID |
     |-------------------------------------|
  1. |      48       TX        Texas     1 |
  2. |      06       CA   California     2 |
  3. |      21       KY     Kentucky     3 |
  4. |      13       GA      Georgia     4 |
  5. |      55       WI    Wisconsin     5 |
     +-------------------------------------+

.

You can get a county shapefile from the Census. Zip codes are trickier to work so you should probably look at the Census' ZIP Code Tabulation Areas (ZCTAs).

[Amanda Kreider]

Hi Robert,

I am creating a map for a project using this shapefile (source), which contains boundaries for a set of "consistent" public-use microdata areas (PUMAs) in the United States. I was hoping to display Alaska and Hawaii next to the contiguous United States, similarly to your example above.

The problem I'm encountering is that the data in my shapefile are already transformed using an Albers projection with standard parallels 29.5 and 45.5. I believe your example requires starting with lat/lon data, and defining different standard parallels for AK/HI vs. the rest of the US. Do you have any idea for how I would go about this?

Thank you so much in advance!
-Amanda