How to add Hawaii and Alaska to spatial polygons in R?

Question

How can I add Hawaii and Alaska to the following code (taken from Josh O'Brien's answer here: Latitude Longitude Coordinates to State Code in R)?

library(sp)
library(maps)
library(maptools)

# The single argument to this function, pointsDF, is a data.frame in which:
#   - column 1 contains the longitude in degrees (negative in the US)
#   - column 2 contains the latitude in degrees

latlong2state <- function(pointsDF) {
    # Prepare SpatialPolygons object with one SpatialPolygon
    # per state (plus DC, minus HI & AK)
    states <- map('state', fill=TRUE, col="transparent", plot=FALSE)
    IDs <- sapply(strsplit(states$names, ":"), function(x) x[1])
    states_sp <- map2SpatialPolygons(states, IDs=IDs,
                     proj4string=CRS("+proj=longlat +datum=wgs84"))

    # Convert pointsDF to a SpatialPoints object 
    pointsSP <- SpatialPoints(pointsDF, 
                    proj4string=CRS("+proj=longlat +datum=wgs84"))

    # Use 'over' to get _indices_ of the Polygons object containing each point 
    indices <- over(pointsSP, states_sp)

    # Return the state names of the Polygons object containing each point
    stateNames <- sapply(states_sp@polygons, function(x) x@ID)
    stateNames[indices]
}

# Test the function using points in Alaska (ak) and Hawaii (hi)

ak <- data.frame(lon = c(-151.0074), lat = c(63.0694))
hi <- data.frame(lon = c(-157.8583), lat = c(21.30694))
nc <- data.frame(lon = c(-77.335), lat = c(34.671))


latlong2state(ak)
latlong2state(hi)
latlong2state(nc)

The latlong2state(ak) and latlong2state(hi) code returns NA, but if the code was modified correctly, Alaska and Hawaii would be returned as results.

Any assistance is appreciated!

Rouen

Answer 1

You need to use a map that has the 50 states, the one you are loading using states <- map('state', fill=TRUE, col="transparent", plot=FALSE) does not have Hawaii and Alaska.

You can for example download the 20m US map from here, and unzip it in your current directory. You should then have a folder called cb_2013_us_state_5m in your R current directory.

I've adapted a bit the code you posted, worked for Hawaii and Alsaka, haven't tried other staes.

library(sp)
library(rgeos)
library(rgdal)

# The single argument to this function, pointsDF, is a data.frame in which:
#   - column 1 contains the longitude in degrees (negative in the US)
#   - column 2 contains the latitude in degrees

latlong2state <- function(pointsDF) {
  states <-readOGR(dsn='cb_2013_us_state_5m',layer='cb_2013_us_state_5m')
  states <- spTransform(states, CRS("+proj=longlat"))

  pointsSP <- SpatialPoints(pointsDF,proj4string=CRS("+proj=longlat"))

  # Use 'over' to get _indices_ of the Polygons object containing each point 
  indices <- over(pointsSP, states)
  indices$NAME
}

# Test the function using points in Alaska (ak) and Hawaii (hi)

ak <- data.frame(lon = c(-151.0074), lat = c(63.0694))
hi <- data.frame(lon = c(-157.8583), lat = c(21.30694))

latlong2state(ak)
latlong2state(hi)

Fibered

NicE - thanks for your help! I had to modify the first reference to states to be states50, and then it worked. I've edited the code above accordingly, but if I'm misunderstanding something, please feel free to revert. – Rouen

This doesnt work for me -- maybe you forget to define states50? – Eulogium

Answer 2

That was based on a data set in the package maps which only included the lower 48. For your task, it needs a shapefile with all the states. Census.gov website is always a good place to find these. I make some changes to the function you posted, so that it works with this new shapefile.

#download a shapefile with ALL states
tmp_dl <- tempfile()
download.file("http://www2.census.gov/geo/tiger/GENZ2013/cb_2013_us_state_20m.zip", tmp_dl)
unzip(tmp_dl, exdir=tempdir())
ST <- readOGR(tempdir(), "cb_2013_us_state_20m")

latlong2state <- function(pointsDF) {
    # Just copied the earlier code with some key changes
    states <- ST

    # Convert pointsDF to a SpatialPoints object 
    # USING THE CRS THAT MATCHES THE SHAPEFILE
    pointsCRS <- "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"
    pointsSP <- SpatialPoints(pointsDF, proj4string=CRS(pointsCRS))

    # Use 'over' to get _indices_ of the Polygons object containing each point 
    indices <- over(pointsSP, states)

    # Return the state names of the Polygons object containing each point
    as.vector(indices$NAME)
}

Let's test it!

ak <- data.frame(lon = c(-151.0074), lat = c(63.0694))
hi <- data.frame(lon = c(-157.8583), lat = c(21.30694))
nc <- data.frame(lon = c(-77.335), lat = c(34.671))

latlong2state(ak)
[1] "Alaska"

latlong2state(hi)
[1] "Hawaii"

latlong2state(nc)
[1] "North Carolina"

Eulogium

J. Winchester - thanks for your assistance. When I run the above code, I receive the following: Error in nchar(projargs) : no method for coercing this S4 class to a vector. If you have any thoughts, please let me know. Thank you. – Rouen

That would seem to be rising from the crs function not recognizing states as the proper class, when I try to extract the projection information from it. I don't know why you get it, as it worked for me. But I modified the answer to work around it. – Eulogium