Search code examples
rr-sp

Points not added to the map plot


I need to add some points to the map using simple points function. The issue is that points don't add to the map. It's simple command, I follow some tutorial where adding points to the map works this way but not in my case. Plot function plots Texas choropleth properly but next line (points) doesn't add points to the map at all:

library(rgdal)
library(rgeos)
library(sp)

companies <- read.csv('geoloc_data_comp.csv', header = T, dec = ',', sep = ';')
states <- readOGR('.', 'states')

plot(states[states@data$stat_name == 'texas',])
points(companies$coords.x1, companies$coords.x2, pch = 21)

Solution

  • First you shoud start to avoid rgeos/rgdal because they will stop being maintains. See : https://github.com/r-spatial/evolution

    sf is replacing them:

    library(sp)
    library(sf)
    library(spData) #used because I wanted US states
    # list of data in spData you have one with US states
    data(package = "spData")
    

    if you want to read shapefile or other GIS format check sf::st_read() (instead of readOGR())

    # one way with sf
    plot(us_states$geometry[us_states$NAME == "Texas"]) 
    # if you want do use the sp way  
    us_sp <- as(us_states, "Spatial") # convert to sp
    plot(us_sp[us_sp@data$NAME == "Texas",])
    

    with sf you have the geometry in one column (see "geometry") instead of having an R S4 with nested lists (see @data and @polygones).

    Before getting some points we need to check in which CRS our data are. If you do not know CRS I like this website : https://ihatecoordinatesystems.com/

    You also have information in the us_states documentation: https://www.rdocumentation.org/packages/spData/versions/2.0.1/topics/us_states

    Then you can use:

    sp::proj4string(us_sp)
    sf::st_crs(us_states)
    # This is EPSG 4269 or NAD83
    

    If you want to use points() they need to be in this coordinates system (I suspect this explain your trouble ie different CRS).

    You didn't provide data points so I produced some:

    library(osmdata)
    #this will just download node matching the key/value place=city
    some_city_in_texas <- osmdata::opq(osmdata::getbb("Texas US"),
                                       nodes_only = TRUE) %>% 
                osmdata::add_osm_feature(key = "place", value = "city") %>% 
                osmdata::osmdata_sf() #keep them in sf format 
                                      # osmdata_sp() also exist
    

    The class osmdata is a bit complicated but here you just need to know that some_city_in_texas$osm_points provide us with points (to test points()). Now we can check their CRS:

    sf::st_crs(some_city_in_texas$osm_points)
    

    As you can see we are in an other CRS so we need to transform it. (you will probably need to do it).

    city_in_texas <- sf::st_transform(some_city_in_texas$osm_points,
                                      4269)
    

    sf use simple feature standard to store localization and points() want two vectors x&y. You should also check that (common cause of error): R use x/y (long/lat) and not lat/long.

    Here we convert city_in_texas to just coords. (if you need to do the reverse, ie converting data frame with X/Y, into an sf object look at sf::st_as_sf())

    coords_city <- sf::st_coordinates(city_in_texas)
    

    Finally this works fine now:

    plot(us_states$geometry[us_states$NAME == "Texas"]) 
    points(coords_city, pch = 21)
    

    Good ressources are https://r-spatial.org/ and https://geocompr.robinlovelace.net/