r - 如何将夏威夷和阿拉斯加添加到 R 中的空间多边形中？

Question

如何将夏威夷和阿拉斯加添加到以下代码中(取自 Josh O'Brien 的回答: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)

latlong2state(ak)和 latlong2state(hi)代码返回 NA，但如果代码修改正确，阿拉斯加和夏威夷将作为结果返回。

任何帮助表示赞赏!

Answer 1

您需要使用具有 50 个州的 map ，您使用 states <- map('state', fill=TRUE, col="transparent", plot=FALSE) 加载的 map 。没有夏威夷和阿拉斯加。

例如，您可以从 here 下载 20m 美国 map 。 , 并将其解压缩到您的当前目录中。然后你应该有一个名为 cb_2013_us_state_5m 的文件夹。在您的 R 当前目录中。

我对您发布的代码进行了一些修改，在夏威夷和阿尔萨卡工作，还没有尝试过其他状态。

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)