https://stat.ethz.ch/pipermail/r-sig-geo/2019-February/027117.html
Reply, then code
You are very close, just needed the expert-knowledge for the incantation required to transform the longlat points (as was already done for the prj.coord, btw).
HTH: https://gist.github.com/mdsumner/afd9228aa1dc4652f8bf193cad4245aa
It does come as a surprise to many that R's plotting is not unit or projection aware (well, it is sometimes, but not usually ... a modern example is ggplot2::geom_sf and the sf package but you need quite a lot of prior-won expertise there too, good luck).
Cheers, Mike.
library(ggplot2)
library(rgdal)
#> Loading required package: sp
#> rgdal: version: 1.3-6, (SVN revision 773)
#> Geospatial Data Abstraction Library extensions to R successfully loaded
#> Loaded GDAL runtime: GDAL 2.3.2, released 2018/09/21
#> Path to GDAL shared files: /usr/share/gdal
#> GDAL binary built with GEOS: TRUE
#> Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
#> Path to PROJ.4 shared files: (autodetected)
#> Linking to sp version: 1.3-1
load(url("https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true"))
PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
NE_countries_rob <- spTransform(NE_countries, CRSobj = PROJ)
NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
NE_box_rob <- spTransform(NE_box, CRSobj = PROJ)
# project long-lat coordinates for graticule label data frames
# (two extra columns with projected XY are created)
prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
lbl.Y.prj <- cbind(prj.coord, lbl.Y)
names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")
prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
lbl.X.prj <- cbind(prj.coord, lbl.X)
names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")
df <-
structure(list(lon = c(2.67569724621467, 17.5766416259819,
28.4126232192772,
23.8147674538232, 29.8917589327105), lat = c(28.1503115976162,
-12.3388787380201, 9.78891068739477, -22.1873831176644, -3.36546931479253
)), class = "data.frame", row.names = c(NA, -5L))
## MDS add projected x,y columns to the existing lon,lat
## - awkward because you can only project a matrix (and there's several other ways ...)
df[c("x", "y")] <- rgdal::project(as.matrix(df[c("lon", "lat")]), PROJ)
m <- ggplot() +
geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
colour="black", fill="gray80", size = 0.25) +
geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black",
fill="transparent", size = 0.25) +
# add graticules projected to Robinson
geom_path(data=NE_graticules_rob, aes(long, lat, group=group),
linetype="dotted", color="grey50", size = 0.25) +
# add graticule labels - latitude and longitude
geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl),
color="grey50", size=2) +
geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl),
color="grey50", size=2) +
coord_fixed(ratio = 1) +
## MDS use our newly projected x,y of the points
geom_point(data=df,
aes(x=x, y=y), colour="Deep Pink",
fill="Pink",pch=21, size=2, alpha=I(0.7)) + theme_void()
#> Regions defined for each Polygons
#> Regions defined for each Polygons
m
Created on 2019-02-18 by the reprex package (v0.2.1)