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

Question

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

5

我该如何在以下代码中添加夏威夷和阿拉斯加（取自Josh O'Brien的答案：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，但如果代码正确修改，则应将阿拉斯加州和夏威夷州作为结果返回。

感谢您的任何帮助！

- Adam Smith

2个回答

2

这是基于包含下48个州的maps数据集的。对于您的任务，需要一个包含所有州的shapefile文件。Census.gov网站通常是找到这些文件的好地方。我对您发布的函数进行了一些更改，以便它可以与这个新的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)
}

让我们来测试一下吧！

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"

- J. Win.

J. Winchester - 感谢您的帮助。当我运行上述代码时，我收到以下错误：Error in nchar(projargs) : no method for coercing this S4 class to a vector。如果您有任何想法，请告诉我。谢谢。 - Adam Smith

似乎是由于crs函数无法识别states作为正确的类而引起的问题，当我尝试从中提取投影信息时。我不知道你为什么会得到它，因为它对我有效。但我修改了答案以解决这个问题。 - J. Win.

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- NicE · Accepted Answer

您需要使用一个包含50个州的地图，您使用的states <- map('state', fill=TRUE, col="transparent", plot=FALSE)加载的地图并不包括夏威夷和阿拉斯加。

例如，您可以从这里下载20m美国地图，并在当前目录中解压缩它。然后，在您的R当前目录中应该有一个名为cb_2013_us_state_5m的文件夹。

我已经修改了您发布的代码，适用于夏威夷和阿拉斯加，但没有尝试其他州。

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)