2

I'm working with GPS data in R from different studies and would like to list the countries that each study covers. I can easily map the data using leaflet but I'm not sure how to extract the countries.

Hope you can help

#' load in the data

mydata <- mydata %>%
  convert(num(long:lat))
str(mydata)

#' map it

 m <- leaflet(mydata) %>%
  addTiles()  %>%
  # setView( lat=-27, lng=170 , zoom=4) %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addCircleMarkers(
    ~ long,
    ~ lat,
    fillOpacity = 0.7,
    color = ~ pal(species),
    radius = 3,
    stroke = FALSE
  ) %>%
  addLegend(
    "bottomright",
    pal = pal,
    values = ~ species,
    title = "Species",
    opacity = 1
  )
m 
3

First you need a map of world polygons.

> library(rworldmap)
> m = getMap()

Then you need your points in a spatial format. I'll make some random points on the globe and create a SpatialPointsDataFrame using the sp package and making sure they have the same projection string as the map:

> set.seed(123); d = data.frame(long=runif(100,-180,180), lat=runif(100,-90,90))
> library(sp)
> coordinates(d) = ~long+lat
> proj4string(d) = proj4string(m)

Now I can extract any of the columns from the map data corresponding to each point with over:

> over(d,m)$NAME
  [1] Jamaica       <NA>          <NA>          <NA>          <NA>         
  [6] <NA>          <NA>          <NA>          Angola        <NA>         
 [11] <NA>          <NA>          Antarctica    <NA>          <NA>         

where <NA> are points in the ocean.

I can add the NAME to the source points:

> d$country = over(d,m)$NAME
> head(d)
             coordinates country
1  (-76.47209, 17.99801) Jamaica
2  (103.7898, -30.09176)    <NA>
3 (-32.76831, -2.049654)    <NA>
4   (137.8863, 81.80529)    <NA>
5  (158.5682, -3.077569)    <NA>
6  (-163.5997, 70.26304)    <NA>

The oceans are big.

Note that: exact world boundaries would be ridiculously complex detailed coastline and border lines. This data set is a simplification. If your points are near the coast they might get mapped to the oceans, or if near a border they might get mapped to the wrong country. Do not use for military planning purposes. Also this is a snapshot of the world at a particular time, countries change boundaries and sovereignties.

If you don't have zillions of these to do and can afford to wait for a web service to give you the answers, use geonames - you will have to get a free geonames user id from geonames.org and then loop over your data:

> library(geonames)
> options(geonamesUsername="myidhere")
> GNcountryCode(17.99801, -76.47209)
$languages
[1] "en-JM"

$distance
[1] "0"

$countryCode
[1] "JM"

$countryName
[1] "Jamaica"
  • That's fantastic, just one more thing. If I run levels(over(d,m)$NAME) it gives me all of the countries in the world. Can I exclude those that have no overlap? – adkane Oct 1 at 16:24
  • I got it. I need to use droplevels. Thanks again – adkane Oct 1 at 16:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.