2013-01-15 3 views
13

ich eine Liste von Breiten- und Längenkoordinaten haben und wollen, um herauszufinden, in welchem ​​Land sie alle wohnen in.Convert Breiten- und Längenkoordinaten zu Ländernamen in R

ich eine Antwort von this question about lat-long to US states modifiziert und haben eine Arbeits Funktion, aber ich laufe auf das Problem, dass die worldHires Karte (aus dem mapdata Paket) schrecklich veraltet ist und viele veraltete Länder wie Jugoslawien und die UdSSR enthält.

Wie würde ich diese Funktion ändern, um ein moderneres Paket zu verwenden, z. B. rworldmap? Ich habe es geschafft mich nur vereiteln so weit ...

library(sp) 
library(maps) 
library(rgeos) 
library(maptools) 

# The single argument to this function, points, is a data.frame in which: 
# - column 1 contains the longitude in degrees 
# - column 2 contains the latitude in degrees 
coords2country = function(points) 
{ 
    # prepare a SpatialPolygons object with one poly per country 
    countries = map('worldHires', fill=TRUE, col="transparent", plot=FALSE) 
    names = sapply(strsplit(countries$names, ":"), function(x) x[1]) 

    # clean up polygons that are out of bounds 
    filter = countries$x < -180 & !is.na(countries$x) 
    countries$x[filter] = -180 

    filter = countries$x > 180 & !is.na(countries$x) 
    countries$x[filter] = 180 

    countriesSP = map2SpatialPolygons(countries, IDs=ids, proj4string=CRS("+proj=longlat +datum=wgs84")) 


    # convert our list of points to a SpatialPoints object 
    pointsSP = SpatialPoints(points, proj4string=CRS("+proj=longlat +datum=wgs84")) 


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


    # Return the state names of the Polygons object containing each point 
    myNames = sapply([email protected], function(x) [email protected]) 
    myNames[indices] 
} 

## 
## this works... but it has obsolete countries in it 
## 

# set up some points to test 
points = data.frame(lon=c(0, 5, 10, 15, 20), lat=c(51.5, 50, 48.5, 47, 44.5)) 

# plot them on a map 
map("worldHires", xlim=c(-10, 30), ylim=c(30, 60)) 
points(points$lon, points$lat, col="red") 

# get a list of country names 
coords2country(points) 
# returns [1] "UK"   "Belgium" "Germany" "Austria" "Yugoslavia" 
# number 5 should probably be in Serbia... 
+0

Die Karten-Paket mit neuen Ländern wird jetzt aktualisiert. Keine UdSSR oder Jugoslawien, etc. – ZacharyST

Antwort

17

Danke für die sorgfältig konstruierte Frage. Es brauchte nur ein paar Zeilenänderungen, um rworldmap (mit aktuellen Ländern) verwenden zu können, siehe unten. Ich bin kein Experte für CRS, aber ich denke nicht, dass die Änderung, die ich an der Projektleitung vornehmen musste, einen Unterschied macht. Andere mögen das gerne kommentieren.

Das ist für mich gearbeitet & gab:

> coords2country(points) 
[1] United Kingdom  Belgium   Germany   Austria   
[5] Republic of Serbia 

die beste All, Andy

library(sp) 
library(rworldmap) 

# The single argument to this function, points, is a data.frame in which: 
# - column 1 contains the longitude in degrees 
# - column 2 contains the latitude in degrees 
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 

    # convert our list of points to a SpatialPoints object 

    # pointsSP = SpatialPoints(points, proj4string=CRS(" +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")) 

    #setting CRS directly to that from rworldmap 
    pointsSP = SpatialPoints(points, proj4string=CRS(proj4string(countriesSP))) 


    # 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 # returns the ISO3 code 
    #indices$continent # returns the continent (6 continent model) 
    #indices$REGION # returns the continent (7 continent model) 
} 
+1

Wenn ich versuche, diese Funktion zu verwenden, bekomme ich 'identicalCRS (x, y) ist nicht wahr ' –

+0

Das sollte durch diese Bearbeitung behoben werden: pointsSP = SpatialPoints (Punkte, proj4string = CRS (proj4string (LänderSP))) – Andy

+0

@Andy könnte dies geändert werden, um einen Kontinent zu extrahieren? und ich habe auch ein paar NA's, als ich das gemacht habe, z. (53.225516, -4.132845, NA) (41.524314, -70.669578, NA) - Weißt du warum? – Ell

8

Sie können meine geonames Paket aus dem http://geonames.org/ Dienst zum Nachschlagen verwenden:

> GNcountryCode(51.5,0) 
$languages 
[1] "en-GB,cy-GB,gd" 

$distance 
[1] "0" 

$countryName 
[1] "United Kingdom of Great Britain and Northern Ireland" 

$countryCode 
[1] "GB" 

> GNcountryCode(44.5,20) 
$languages 
[1] "sr,hu,bs,rom" 

$distance 
[1] "0" 

$countryName 
[1] "Serbia" 

$countryCode 
[1] "RS" 

es von r-Schmiede, weil ich bin nicht sicher, dass ich es zu CRAN zu lösen gestört:

https://r-forge.r-project.org/projects/geonames/

Ja, es ist abhängig von einem externen Dienst, aber zumindest weiß, was Happe Ned zum Kommunismus ... :)

+0

Es scheint, dass die API einen Benutzernamen benötigt, der in jeden Aufruf eingebaut wird, aber die 'GNcountryCode' Funktion erlaubt keinen Benutzernamen. Werden Sie die API-Aufrufe anpassen? @ Spacedman –

Verwandte Themen