[R-sig-Geo] Plot coordinates on world map with Robinson CRS - ggplot2

Michael Sumner md@umner @end|ng |rom gm@||@com
Mon Feb 18 11:00:59 CET 2019


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.

On Mon, 18 Feb 2019 at 20:01 Miluji Sb <milujisb using gmail.com> wrote:

> Dear all,
>
> I am trying to plot coordinates on a world map with Robinson CRS. While the
> world map is generated without any issues, when I try to plot the points -
> I only get a single point.
>
> The code I am using and the coordinates data is below. What am I doing
> wrong? Any help/suggestions will be highly appreciated.
>
> library(data.table)
> library(foreign)
> library(readstata13)
> library(rgdal)
> library(maptools)
> library(ggplot2)
> library(dplyr)
>
> 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")
>
> m <- ggplot() +
>   # add Natural Earth countries projected to Robinson, give black border
> and fill with gray
>   geom_polygon(data=NE_countries_rob, aes(long,lat, group=group),
> colour="black", fill="gray80", size = 0.25) +
>   # Note: "Regions defined for each Polygons" warning has to do with
> fortify transformation. Might get deprecated in future!
>   # alternatively, use use map_data(NE_countries) to transform to data
> frame and then use project() to change to desired projection.
>   # add Natural Earth box projected to Robinson
>   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) +
>   # the default, ratio = 1 in coord_fixed ensures that one unit on the
> x-axis is the same length as one unit on the y-axis
>   coord_fixed(ratio = 1) +
>   geom_point(data=df,
>              aes(x=lon, y=lat), colour="Deep Pink",
>              fill="Pink",pch=21, size=2, alpha=I(0.7))
>   # remove the background and default gridlines
>   theme_void()
>
> ## coordinates dataframe
> dput(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))
>
> Sincerely,
>
> Milu
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list