You need to make a spatial points object and then reproject __that__: library(sp) library(maptools) data(wrld_simpl) wrld_simpl.trans <- spTransform(wrld_simpl, CRS = ("+proj=moll +lon_0=0 +x_0=0 +y_0=0")) plot(wrld_simpl.trans, col = 'lightgray')) latitudes <- c(-30, 0, 20, 50) longitudes <- c(-20, 10, 50, 80) Make an SPDF: pts <- data.frame(latitudes, longitudes) coordinates(pts) <- ~longitudes+latitudes proj4string(pts) <- CRS(paste("+init=epsg:4326")) Project and plot it: pts.trans <- spTransform(pts, CRS = ("+proj=moll +lon_0=0 +x_0=0 +y_0=0")) points(pts.trans, cex = 1.5, col = 'red', pch = 19) [![enter image description here][1]][1] # Compare to unprojected plot(wrld_simpl, col = 'lightgray') points(pts, cex = 1.5, col = 'red', pch = 19) [![enter image description here][2]][2] EDIT: This is actually a little more fun with extreme points: longitudes <- seq(-150, 150, length.out = 5) latitudes <- seq(-90, 90, length.out = 5) [![enter image description here][3]][3] [1]: https://i.sstatic.net/soNOI.png [2]: https://i.sstatic.net/NGqsL.png [3]: https://i.sstatic.net/hn7ju.png