-1

我正在使用一个包含超过 270,000 个观察值的三变量数据集。这三个变量是观测值、纬度和经度。在以前的相关帖子中，我设法获得有关如何允许反向地理编码功能跳过纬度和经度缺失值的观察的帮助：使用反向地理编码功能时如何处理缺失值？

可重现的例子：

Data <- data.frame(
  Observation = 1:5,
  Longitude = c(116.3880005, 53, -97, NA, NA), 
  Latitude = c(39.92890167, 32, 32, NA, NA))

以下代码有效。但是，它会为每个国家/地区生成因子指数，而不是我希望获得的 ISO3。

library(sp)
library(rworldmap)

coords2country = function(points)
{  
  countriesSP <- getMap(resolution='low')
  #countriesSP <- getMap(resolution='high') #you could use high res map from rworldxtra if you were concerned about detail
  # new changes in worldmap means you have to use this new CRS (bogdan):
  pointsSP = SpatialPoints(points, proj4string=CRS(" +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))

  # use 'over' to get indices of the Polygons object containing each point 
  indices = over(pointsSP, countriesSP)

  # return the ADMIN names of each country
  indices$ADMIN  
  indices$ISO3 #would return the ISO3 code
}

#The function below fixes the NA error
coords2country_NAsafe <- function(points)
{
  bad <- with(points, is.na(lon) | is.na(lat))
  result <- character(length(bad))
  result[!bad] <- coords2country(points[!bad,])
  result
}

以下生成因子索引（不是 ISO3 代码）：

coords2country_NAsafe(points)

我想知道我可以修改上面的代码以输出 ISO3 代码，而不是它们的因子索引。

4

2 回答 2

2

This all looks good to me:

> pts=read.table("points.csv",head=TRUE,sep=",")
> pts=subset(pts,!is.na(lon)) # just take the good ones to save NA faffing
> coordinates(pts)=~lon+lat
> first = pts[1:100,]         # take first 100 for starters
> cc = coords2country(first)
> plot(first,pch=".")
> text(coordinates(first),label=cc)

enter image description here

All the countries in the right place...

于 2013-08-21T12:38:23.470 回答
2

我认为您的主要问题是您输出的是每个 ISO3 代码的因子索引，而不是 ISO3 代码本身。因此，中国有 42 个，因为中国是地图上的第 42 个国家。下面的 as.character() 对其进行了排序。

因此，对您的 & Barry 的代码进行少量编辑会给出下面的代码，我认为它可以满足您的需求。

在最后 4 行中将 'first' 替换为 'pts' 以运行整个数据集。

coords2country = function(points)
{  
  library(rworldmap)
  countriesSP <- getMap(resolution='low')

  #I assume that points have same projection as the map
  pointsSP = SpatialPoints(points, proj4string=CRS(proj4string(countriesSP)))  

  # use 'over' to get indices of the Polygons object containing each point 
  indices = over(pointsSP, countriesSP)

  #as.character(indices$ADMIN) # return the ADMIN names of each country  
  as.character(indices$ISO3) # return the ISO3 code
  #as.character(indices$ISO_N3) # return the ISO numeric code
}


library(sp)
pts=read.table("points.csv",head=TRUE,sep=",")
pts=subset(pts,!is.na(lon)) # just take the good ones to save NA faffing
coordinates(pts)=~lon+lat
first = pts[1:100,]         # take first 100 for starters
cc = coords2country(first)
plot(first,pch=".")
text(coordinates(first),label=cc)

firstPlusISO3 <- cbind(coordinates(first),cc)
于 2013-08-21T21:54:34.363 回答