Как получить код ISO3, а не индекс фактора, из моей функции обратного геокода?

Question

Я работаю с набором данных из трех переменных с более чем 270000 наблюдений. Три переменные - это значение наблюдения, широта и долгота. В предыдущем связанном посте мне удалось получить справку о том, как разрешить функции обратного геокода пропускать наблюдения с отсутствующими значениями широты и долготы: как обрабатывать пропущенные значения при использовании функции обратного геокода?

Воспроизводимый пример:

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, а не их факторные индексы.

Answer 1

Я думаю, что ваша основная проблема заключалась в том, что вы выводили индекс фактора для каждого кода ISO3 вместо самих кодов ISO3. Таким образом, у вас было 42 для Китая, потому что Китай является 42-й страной на карте. As.character() ниже сортирует это.

Таким образом, внесение небольших изменений в код вашего & Барри дает код, который, как я думаю, делает то, что вы хотите.

Замените "first" на "pts" в последних 4 строках, чтобы запустить весь ваш набор данных.

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)  

Answer 2

Это все выглядит хорошо для меня:

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

Все страны в нужном месте...