2013-01-15 16 views
13

Mam listę szerokości i długości geograficznej, i chcą się dowiedzieć, w jakim kraju wszyscy mieszkają w.Konwersja szerokość i długość geograficzną nazwę państwa na badania

I zmodyfikowana odpowiedź od this question about lat-long to US states i mają pracę funkcja, ale napotkam problem, że mapa worldHires (z pakietu mapdata) jest ohydnie nieaktualna i zawiera wiele przestarzałych krajów, takich jak Jugosławia i ZSRR.

Jak zmodyfikować tę funkcję, aby użyć bardziej nowoczesnego pakietu, takiego jak rworldmap? Mam tylko udało się udaremnić siebie tak daleko ...

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

Pakiet map jest obecnie aktualizowany o nowe kraje. Nie ma już ZSRR ani Jugosławii, itp. – ZacharyST

Odpowiedz

17

Dzięki za starannie skonstruowane pytanie. Wymagało tylko kilku zmian linii, aby móc korzystać z rworldmap (zawierającego aktualne kraje), patrz poniżej. Nie jestem ekspertem od CRS, ale nie sądzę, żeby zmiana, jaką musiałem wprowadzić w proj4stringach, miała znaczenie. Inni mogą chcieć to skomentować.

ten pracował dla mnie & dał:

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

Wszystko co najlepsze, 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

Kiedy próbuję użyć tej funkcji, otrzymuję 'identyczneCRS (x, y) nie jest PRAWDZIWE' –

+0

To powinno być naprawione przez tę edycję: pointsSP = SpatialPoints (points, proj4string = CRS (proj4string (countriesSP))) – Andy

+0

@Andy można to zmodyfikować, aby wyodrębnić kontynent? i dostałem też trochę NA, kiedy to uruchomiłem, np. (53.225516, -4.132845, NA) (41,524314, -70,669578, NA) - czy wiesz dlaczego? – Ell

8

Można użyć paczkę geonames do wyszukiwania z serwisu http://geonames.org/:

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

Pobierz go z R-kuźni, bo nie jestem pewien, że jedno, aby zwolnić go do CRAN:

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

Tak, to zależy od służby zewnętrznej, ale przynajmniej wie co Happe ned do komunizmu ... :)

+0

Wygląda na to, że API wymaga nazwy użytkownika wbudowanej w każde połączenie, ale funkcja 'GNcountryCode' nie zezwala na nazwę użytkownika. Czy będziesz dostosowywać wywołania API? @Spacedman –

Powiązane problemy