US Maps using Winston Chang's R Graphic Cookbook (pages 317-318)

967 views
Skip to first unread message

pradip...@samhsa.hhs.gov

unread,
Jan 19, 2013, 8:25:36 PM1/19/13
to ggp...@googlegroups.com

Hello,

My goal is to draw a US map similar to the example provided in Winston Chang's R Graphic Cookbook (pages 317-318) except that I wanted to add a title and add Alaska and Hawaii to the map.
The reproducible example is appended below .

There are 2 issues:

1) I got the following error message after the theme_clean() function:

Error in seq_len(nrow(data) - 1) : 
  argument must be coercible to non-negative integer

2) I could not figure out how to add Alaska and Hawaii to the map.

I would appreciate if someone could help me resolve the issuea.

Thanks in advance.

Pradip Muhuri

The reproducible example begins here

#################
setwd ("E:/CSAP_P1")
library(maps)
library(plyr)
loadedNamespaces()

#read estimates of risk factors for substances use (ages 12-17) by State obtained from SUDAAN output
p1_st_data <-read.table (text="
Alabama,49.60,    1.37
Alaska,55.00,    1.41
Arizona,52.50,    1.56
Arkansas,50.50,    1.22
California,    51.10,    0.65
Colorado,55.10,    1.26
Connecticut,56.30,    1.28
Delaware,    53.60,    1.30
District of Columbia,  53.50,    1.22
Florida,  52.70,    0.67
Georgia,    52.50,    1.15
Hawaii,    49.40,    1.33
Idaho,    48.30,    1.23
Illinois,    52.70,    0.63
Indiana,    49.60,    1.16
Iowa,    46.30,    1.37
Kansas,    44.30,    1.43
Kentucky,    52.90,    1.37
Louisiana,    49.70,    1.23
Maine,    55.60,    1.44
Maryland,    53.90,    1.46
Massachusetts,    55.40,    1.41
Michigan,    52.40,    0.62
Minnesota,    51.50,    1.20
Mississippi,    43.20,    1.14
Missouri,    48.70,    1.20
Montana,    56.40,    1.16
Nebraska,    45.70,    1.51
Nevada,    54.20,    1.17
New Hampshire,    56.10,    1.30
New Jersey,    53.20,    1.45
New Mexico,    57.60,    1.34
New York,    53.70,    0.67
North Carolina,    52.20,    1.26
North Dakota,    48.60,    1.34
Ohio,    50.90,    0.61
Oklahoma,    47.20,    1.42
Oregon,    54.00,    1.35
Pennsylvania,    53.00,    0.63
Rhode Island,    57.20,    1.20
South Carolina,    50.50,    1.21
South Dakota,    43.40,    1.30
Tennessee,    48.90,    1.35
Texas,    48.70,    0.62
Utah,    42.00,    1.49
Vermont,    58.70,    1.24
Virginia,    51.80,    1.18
Washington,    53.50,    1.39
West Virginia,    52.80,    1.07
Wisconsin,    49.90,    1.50
Wyoming,    49.20,    1.29",
 sep=  "," , col.names = c("state" ,   "Obt_mrj_p" ,  "Obt_mrj_se" ),
 colClasses = c( "character" ,  "numeric" , "numeric" )                     
)   

#create a new variable, region
names(p1_st_data) <- tolower (names(p1_st_data))
p1_st_data$region <- tolower(p1_st_data$state)

p1_st_data$ob_mrj_cat <- cut (p1_st_data$obt_mrj_p, quantile (p1_st_data$obt_mrj_p, (0:5/5)),
                                                  include.lowest=TRUE)


#load maps data library
library(maps)
states_map <- map_data("state")

head (states_map)

#merge US map data with NSDUH estimates by region (state)
drug_map <- merge (states_map, p1_st_data,
                   by.states_map = "region", by.p1_st_data = "state")

head (drug_map)


#sort by group, then order
drug_map <- arrange (drug_map, group, order)


theme_clean <- function (base_size=12) {
  require(gridBase)
  theme_grey(base_size) %+replace%
   
    theme(
    axis.title = element_blank (),
    axis.text = element_blank (),
    panel.background = element_blank (),
    panel.grid = element_blank (),
    axis.ticks.length = unit (0,"cm"),
    axis.ticks.margin = unit (0,"cm"),
    panel.margin = unit (0,"lines"),
    plot.margin = unit(c(0,0,0,0),"lines"),
    complete =TRUE
  )
   
}

#generate  duscrete color palete with 5 values
pal <- colorRampPalette(c("#559999", "grey80", "#BB650B")) (5)

# US map on % of those who said - easy to get marijuana (Obt_mrj_p)
 ggplot(drug_map, aes(map_id=state, fill = ob_mrj_cat)) +
                 geom_map (map=states_map, colour="black")+
                  scale_fill_manual (values=pal)+
                  expand_limits (x=states_map$long, y=states_map$lat)+
                  coord_map ("polyconic")+
                  labs (fill= "Pecentages")+
  ggtitle ("Fairly Easy/Very Easy to Obtain Marijuana
           \nas Reported by Youths (Ages 12 to 17)
           \nPercentages, Annual Averages 2002-2008 NSDUH  ")+
            
theme_clean()

ggsave (file="US_map_mrj_rev.png", width=11, height=7)



Winston Chang

unread,
Jan 19, 2013, 9:40:45 PM1/19/13
to pradip...@samhsa.hhs.gov, ggplot2
Hi Pradip -

The issue is that you need to use map_id=region instead of map_id=state, because in your two data sets, that's the name of the column that identifies the state. (You do have a 'state' column in drug_map, but not in states_map.)

Also, some minor things with theme_clean:
- Use require(grid), not require(gridBase)
- Need to use axis.ticks.margin=unit(0.01,"cm") instead of unit(0,"cm"). As noted in the text, this is a workaround for a bug that's present in R<=2.15.2, but will be fixed in the next released version of R.


For information about adding Alaska and Hawaii, please see:

The code (starting at the call to theme_clean) should be:

theme_clean <- function (base_size=12) {
  require(grid) 
  theme_grey(base_size) %+replace%
    theme(
      axis.title = element_blank (),
      axis.text = element_blank (),
      panel.background = element_blank (),
      panel.grid = element_blank (),
      axis.ticks.length = unit (0,"cm"),
      axis.ticks.margin = unit (0.01,"cm"),
      panel.margin = unit (0,"lines"),
      plot.margin = unit(c(0,0,0,0),"lines"),
      complete = TRUE
    )    
}

#generate  duscrete color palete with 5 values
pal <- colorRampPalette(c("#559999", "grey80", "#BB650B")) (5)

# US map on % of those who said - easy to get marijuana (Obt_mrj_p)
ggplot(drug_map, aes(map_id=region, fill = ob_mrj_cat)) + 
  geom_map (map=states_map, colour="black")+
  scale_fill_manual (values=pal)+
  expand_limits (x=states_map$long, y=states_map$lat)+
  coord_map ("polyconic")+
  labs (fill= "Pecentages")+
  ggtitle ("Fairly Easy/Very Easy to Obtain Marijuana
           \nas Reported by Youths (Ages 12 to 17)
           \nPercentages, Annual Averages 2002-2008 NSDUH  ")+
  theme_clean()


=========================================================
As for adding Alaska and Hawaii, I've adapted this answer from stackoverflow.

At the end, it just makes a map with all the states, but since the data set doesn't have state names in them, you'll have to do some modifications to add the names to the data. Then you should be able to combine the map data with your data of interest.

Note that the resulting map already has a projection applied, so you don't need to use coord_map("polyconic") with it.


require(maptools)
require(rgdal)

fixup <- function(usa,alaskaFix,hawaiiFix){

  alaska=usa[usa$STATE_NAME=="Alaska",]
  alaska = fix1(alaska,alaskaFix)
  proj4string(alaska) <- proj4string(usa)

  hawaii = usa[usa$STATE_NAME=="Hawaii",]
  hawaii = fix1(hawaii,hawaiiFix)
  proj4string(hawaii) <- proj4string(usa)

  usa = usa[! usa$STATE_NAME %in% c("Alaska","Hawaii"),]
  usa = rbind(usa,alaska,hawaii)

  return(usa)
}

fix1 <- function(object,params){
  r=params[1];scale=params[2];shift=params[3:4]
  object = elide(object,rotate=r)
  size = max(apply(bbox(object),1,diff))/scale
  object = elide(object,scale=size)
  object = elide(object,shift=shift)
  object
}

# Use downloader library to reliably download https across platforms
library(downloader)
         "usmapdata.zip")
# File originally from:

# Unzip. This zip file happens to have extraneous files with __MACOSX in path,
# so ignore them as we unzip.
files <- as.character(unzip("usmapdata.zip", list=TRUE)$Name)
files <- files[!grepl("__MACOSX", files)]
unzip("usmapdata.zip", files = files)

# Read in shapefile using rgdal
us <- readOGR(dsn = "states_21basic", layer = "states")

# Transform to equal area and run fixup
usAEA = spTransform(us,CRS("+init=epsg:2163"))
usfix = fixup(usAEA,c(-35,1.5,-2800000,-2600000),c(-35,1,6800000,-1600000))

# Convert to data frame
library(ggplot2)
usfix_df <- fortify(usfix)
head(usfix_df)
#       long      lat order  hole piece group id
# 1 -1632777 590440.6     1 FALSE     1   1.1  1
# 2 -1636961 592035.6     2 FALSE     1   1.1  1
# 3 -1639547 581112.1     3 FALSE     1   1.1  1
# 4 -1635144 571825.2     4 FALSE     1   1.1  1
# 5 -1643611 582947.5     5 FALSE     1   1.1  1
# 6 -1643111 591750.0     6 FALSE     1   1.1  1

# Plot
ggplot(usfix_df, aes(x=long, y=lat, group=group)) + geom_polygon()


I hope this helps!

-Winston








--
You received this message because you are subscribed to the ggplot2 mailing list.
Please provide a reproducible example: https://github.com/hadley/devtools/wiki/Reproducibility
 
To post: email ggp...@googlegroups.com
To unsubscribe: email ggplot2+u...@googlegroups.com
More options: http://groups.google.com/group/ggplot2

Muhuri, Pradip (SAMHSA/CBHSQ)

unread,
Jan 21, 2013, 9:50:28 AM1/21/13
to Winston Chang, ggplot2, Muhuri, Pradip (SAMHSA/CBHSQ)

Hi Winston,

 

Thank you so much for your help with the revision of my code.  The revised reproducible example is attached.

 

When running the revised code, I am getting the following message. Why am I getting this error message?  Any thoughts?

 

 

> #load maps data library

> library(maps)

> states_map <- map_data("region")

Error in get(dbname) : object 'regionMapEnv' not found

In addition: Warning message:

In data(list = dbname) : data set ‘regionMapEnv’ not found

 

 

http://finzi.psych.upenn.edu/R/library/ggplot2/html/map_data.html

 

 

Congratulations on writing the R Graphics Cookbook! I find the examples in the book very helpful.

 

 

regards,

 

Pradip

 

 

Pradip K. Muhuri, PhD

Statistician

Substance Abuse & Mental Health Services Administration

The Center for Behavioral Health Statistics and Quality

Division of Population Surveys

1 Choke Cherry Road, Room 2-1071

Rockville, MD 20857

 

Tel: 240-276-1070

Fax: 240-276-1260

e-mail: Pradip...@samhsa.hhs.gov

 

The Center for Behavioral Health Statistics and Quality your feedback.  Please click on the following link to complete a brief customer survey:   http://cbhsqsurvey.samhsa.gov

drug_map_rev.R

Rainer Hurling

unread,
Jan 21, 2013, 3:46:26 PM1/21/13
to Muhuri, Pradip (SAMHSA/CBHSQ), Winston Chang, ggplot2
On 21.01.2013 15:50 (UTC+2), Muhuri, Pradip (SAMHSA/CBHSQ) wrote:
> Hi Winston,
>
>
>
> Thank you so much for your help with the revision of my code. The
> revised reproducible example is attached.
>
>
>
> When running the revised code, I am getting the following message. Why
> am I getting this error message? Any thoughts?
>
>
>
>
>
>> #load maps data library
>
>> library(maps)
>
>> states_map <- map_data("region")

Hi Pradip,

your attached file again differs from your inlined example by having
map_data("region") instead of map_data("state"), as Winston proposed
before. So the error is exactly, what its message tells us: no Map
Environment for "region".

HTH,
Rainer

> Error in get(dbname) : object 'regionMapEnv' not found
>
> In addition: Warning message:
>
> In data(list = dbname) : data set �regionMapEnv� not found
>
>
>
>
>
> http://finzi.psych.upenn.edu/R/library/ggplot2/html/map_data.html
>
>
>
>
>
> Congratulations on writing the R Graphics Cookbook! I find the examples
> in the book very helpful.
>
>
>
>
>
> regards,
>
>
>
> Pradip
>
>
>
>
>
> Pradip K. Muhuri, PhD
>
> Statistician
>
> Substance Abuse & Mental Health Services Administration
>
> The Center for Behavioral Health Statistics and Quality
>
> Division of Population Surveys
>
> 1 Choke Cherry Road, Room 2-1071
>
> Rockville, MD 20857
>
>
>
> Tel: 240-276-1070
>
> Fax: 240-276-1260
>
> e-mail: Pradip...@samhsa.hhs.gov <mailto:Pradip...@samhsa.hhs.gov>
>
>
>
> The Center for Behavioral Health Statistics and Quality your feedback.
> Please click on the following link to complete a brief customer
> survey: http://cbhsqsurvey.samhsa.gov <http://cbhsqsurvey.samhsa.gov/>
>
>
>
> *From:*Winston Chang [mailto:winsto...@gmail.com]
> *Sent:* Saturday, January 19, 2013 9:41 PM
> *To:* Muhuri, Pradip (SAMHSA/CBHSQ)
> *Cc:* ggplot2
> *Subject:* Re: US Maps using Winston Chang's R Graphic Cookbook (pages
> <mailto:pradip...@samhsa.hhs.gov> <pradip...@samhsa.hhs.gov

Muhuri, Pradip (SAMHSA/CBHSQ)

unread,
Jan 22, 2013, 5:51:12 PM1/22/13
to Rainer Hurling, Winston Chang, ggplot2
Hi Winston and Rainer,

Thank you so much for your help. The revised code worked fine.

Regards,

Pradip


Pradip K. Muhuri, PhD
Statistician
Substance Abuse & Mental Health Services Administration
The Center for Behavioral Health Statistics and Quality
Division of Population Surveys
1 Choke Cherry Road, Room 2-1071
Rockville, MD 20857
 
Tel: 240-276-1070
Fax: 240-276-1260
e-mail: Pradip...@samhsa.hhs.gov
 
The Center for Behavioral Health Statistics and Quality your feedback.  Please click on the following link to complete a brief customer survey:  http://cbhsqsurvey.samhsa.gov

-----Original Message-----
From: Rainer Hurling [mailto:rhu...@gwdg.de]
Sent: Monday, January 21, 2013 3:46 PM
To: Muhuri, Pradip (SAMHSA/CBHSQ)
Cc: 'Winston Chang'; ggplot2
Subject: Re: US Maps using Winston Chang's R Graphic Cookbook (pages 317-318)

On 21.01.2013 15:50 (UTC+2), Muhuri, Pradip (SAMHSA/CBHSQ) wrote:
> Hi Winston,
>
>
>
> Thank you so much for your help with the revision of my code. The
> revised reproducible example is attached.
>
>
>
> When running the revised code, I am getting the following message. Why
> am I getting this error message? Any thoughts?
>
>
>
>
>
>> #load maps data library
>
>> library(maps)
>
>> states_map <- map_data("region")

Hi Pradip,

your attached file again differs from your inlined example by having
map_data("region") instead of map_data("state"), as Winston proposed
before. So the error is exactly, what its message tells us: no Map
Environment for "region".

HTH,
Rainer

> Error in get(dbname) : object 'regionMapEnv' not found
>
> In addition: Warning message:
>
> In data(list = dbname) : data set 'regionMapEnv' not found
US_map_mrj_rev.png
Reply all
Reply to author
Forward
0 new messages