Displaying Alaska and Hawaii counties next to the US using spmap

[Marina Bosque]

Hi everyone,

I am trying to draw a US map with all the counties using the package spmap. The problem is that when Stata plots the map and I include Alaska and Hawaii, the map appears distorted. The data used can be downloaded from here: https://www.census.gov/geo/maps-data..._counties.html.

The code I am using is the following one:

    cd "../cb_2016_us_county_500k"
    ssc install spmap
    ssc install shp2dta

    shp2dta using cb_2016_us_county_500k, database(usdb) coordinates(uscoord) genid(id)

    ****
    use "usdb.dta", clear
    
    destring STATEFP, gen(STATEFP_new)
    drop STATEFP
    rename STATEFP_new STATEFP
    
    rename STATEFP state_fips 
    rename NAME county_str

    replace county_str = subinstr(county_str, "Baltimore", "Baltimore City", .) if LSAD=="25"
    replace county_str = subinstr(county_str, "St. Louis", "St. Louis City", .) if LSAD=="25"
    replace county_str = subinstr(county_str, "Fairfax", "Fairfax City", .) if LSAD=="25"
    replace county_str = subinstr(county_str, "Franklin", "Franklin City", .) if LSAD=="25"
    replace county_str = subinstr(county_str, "Richmond", "Richmond City", .) if LSAD=="25"
    replace county_str = subinstr(county_str, "Roanoke", "Roanoke City", .) if LSAD=="25"
    replace county_str = subinstr(county_str, "Doña Ana", "Dona Ana", .)
    // Rename the variables to merge
    
    ** Merge datasets
    
    merge 1:1 state_fips county_str using "../Crosswalk_county_raw_temp.dta"
    // Merge with the dataset that has the data we want to plot.
    
    drop if _merge!=3    
    drop _merge
    drop if state_fips==72 // Dropping Puerto Rico
    
    ** Draw the graph
    
    spmap color_num using "uscoord.dta", id(id) ///
    fcolor(RdYlBu) ndfcolor(white) osize(thin) clmethod(custom) clbreaks(0 1 2 3 4 5) ///
    legtitle(Belt Definitions) legorder(lohi) legend(label(2 "Rust Belt") label(3 "Sun Belt") ///
    label(4 "East Coast") label(5 "West Coast") label(6 "Flyover counties") position(2)) ///
    graphregion(fcolor(white) lcolor(white) ilcolor(white))  ///
    title({bf:Potential Belt Definitions}) 

With this code, Alaska and Hawaii distort the map. I have been trying to modify the uscoord.dta file in the way it was shown in previous posts (for a state level map), but I don't know if I am doing it well or if I should do it in another way. The error Stata gave me is:

latitude _Y must be between -90 and 90

. The following part of code would go after the command

shp2dta using cb_2016_us_county_500k, database(usdb) coordinates(uscoord) genid(id)

    ssc install geo2xy
    net get geo2xy

    use "../cb_2016_us_county_500k/uscoord.dta", clear

    * Alaska
    drop if _X > 0 & !mi(_X) & !inlist(_ID, 28,29,30,463,464,465,466,467,468,774,900,953,1054,1076,1158,1278,1369,1370,1371,2135,2276,2277,2278,2279,2280,2281,2527,2556,3073)  // land in the Eastern hemisphere
    geo2xy _Y _X if !inlist(_ID, 28,29,30,463,464,465,466,467,468,774,900,953,1054,1076,1158,1278,1369,1370,1371,2135,2276,2277,2278,2279,2280,2281,2527,2556,3073), replace proj(albers)
    replace _Y = _Y / 3 + 700000 if !inlist(_ID, 28,29,30,463,464,465,466,467,468,774,900,953,1054,1076,1158,1278,1369,1370,1371,2135,2276,2277,2278,2279,2280,2281,2527,2556,3073)
    replace _X = _X / 3 - 1700000 if !inlist(_ID, 28,29,30,463,464,465,466,467,468,774,900,953,1054,1076,1158,1278,1369,1370,1371,2135,2276,2277,2278,2279,2280,2281,2527,2556,3073)

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

    * contiguous US
    geo2xy _Y _X if !inlist(_ID,98,359,360,361,2389,28,29,30,463,464,465,466,467,468,774,900,953,1054,1076,1158,1278,1369,1370,1371,2135,2276,2277,2278,2279,2280,2281,2527,2556,3073), replace proj(albers)

The IDs are the ones that correspond to the counties in Hawaii and Alaska according to the usdb.dta file.

[Robert Picard]

It's a lot to ask to plot a map of every US county!

Here's a scoop of an upcoming update to geo2xy that includes a "Fun with Maps" help page where I show how to create a composite map of the US (at the state level).

The following starts from the 500K shapefile from the Census (link in #1). You may want to play with the scale and offsets for Alaska, Hawaii, and Puerto Rico.

shp2dta using cb_2016_us_county_500k/cb_2016_us_county_500k, ///
    database(usdb) coordinates(uscoord) genid(id) replace
    
* merge state identifier with coordinates
use "uscoord.dta", clear
gen long order = _n
tempfile f
save "`f'"
use "usdb.dta", clear
keep STATEFP id
rename id _ID
merge 1:m _ID using "`f'", assert(match) nogen
sort _ID order

* drop US territories
* 60 == American Samoa
* 66 == Guam
* 69 == Commonwealth of the Northern Mariana Islands
* 78 == U.S. Virgin Islands
drop if inlist(STATEFP, "60", "66", "69", "78")

gen byte Alaska = STATEFP == "02"
gen byte Hawaii = STATEFP == "15"
drop if Hawaii & _X < -160 // small islands to the west
gen byte PuertoRico = STATEFP == "72"
gen byte us48 = !(Alaska | Hawaii | PuertoRico)

* flip longitudes to reconnect Hawaii and Alaska
replace _X = cond(_X > 0, _X - 180, _X + 180) if Alaska | Hawaii
    
* Alaska - USGS recommends standard parallels of 55 and 65 north
sum _X if Alaska
local midlon = (r(min) + r(max)) / 2
geo2xy _Y _X if Alaska, replace ///
        proj(albers, 6378137 298.257223563 55 65 0 `midlon')
replace _Y = _Y / 3 + 800000 if Alaska
replace _X = _X / 3 - 1700000 if Alaska

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

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

save "xy_coor.dta", replace
    
use "usdb.dta", clear
spmap using "xy_coor.dta", id(id)

and here's the map produced

[Marina Bosque]

Hi Robert,

Many thanks for your help, it perfectly works!