Search code examples
rmapsgisr-maptoolsr-sp

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


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!


Solution

  • 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)