Plot coordinates on world map with Robinson CRS - ggplot2

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Plot coordinates on world map with Robinson CRS - ggplot2

Miluji Sb
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Roger Bivand
Administrator
On Mon, 18 Feb 2019, Miluji Sb 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.

Well, the point is correct, since you didn't project df, did you? This
seems to be typical of the arcane nature of ggplot invocations. Did you
try tmap? Almost all the loaded&attached packages are superfluous too.

Roger

>
> 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Michael Sumner-2
In reply to this post by Miluji Sb
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 <[hidden email]> 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
> [hidden email]
> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Plot coordinates on world map with Robinson CRS - ggplot2

Miluji Sb
Thank you very much. This works perfectly. Appreciate it.

Sincerely,

Milu

On Mon, Feb 18, 2019 at 11:01 AM Michael Sumner <[hidden email]> wrote:

> 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 <[hidden email]> 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
>> [hidden email]
>> 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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo