-
-
Save valentinitnelav/c7598fcfc8e53658f66feea9d3bafb40 to your computer and use it in GitHub Desktop.
# ====================================================================================== | |
# Pacific centered world map with ggplot | |
# Enhanced aspect with graticules and labels | |
# The central/prime meridian can be shifted with any positive value towards west | |
# Can use any project of known PROJ.4 string | |
# ====================================================================================== | |
# ~~~~~~~~~~~ Load needed libraries ~~~~~~~~~~~ # | |
library(data.table) | |
library(ggplot2) | |
library(rgdal) | |
library(rgeos) | |
library(maps) | |
library(maptools) | |
library(raster) | |
# ~~~~~~~~~~~ Download shapefile from www.naturalearthdata.com ~~~~~~~~~~~ # | |
# Download countries data | |
download.file(url = "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip", | |
destfile = "ne_110m_admin_0_countries.zip") | |
# unzip the shapefile in the directory mentioned with "exdir" argument | |
unzip(zipfile="ne_110m_admin_0_countries.zip", exdir = "ne_110m_admin_0_countries") | |
# delete the zip file | |
file.remove("ne_110m_admin_0_countries.zip") | |
# read the shapefile with readOGR from rgdal package | |
NE_countries <- readOGR(dsn = "ne_110m_admin_0_countries", layer = "ne_110m_admin_0_countries") | |
class(NE_countries) # is a SpatialPolygonsDataFrame object | |
# ~~~~~~~~~~~ Split world map by "split line" ~~~~~~~~~~~ # | |
# inspired from: | |
# https://stat.ethz.ch/pipermail/r-sig-geo/2015-July/023168.html | |
# shift central/prime meridian towards west – positive values only | |
shift <- 180+30 | |
# create "split line" to split country polygons | |
WGS84 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") | |
split.line = SpatialLines(list(Lines(list(Line(cbind(180-shift,c(-90,90)))), ID="line")), | |
proj4string=WGS84) | |
# NOTE - in case of TopologyException' errors when intersecting line with country polygons, | |
# apply the gBuffer solution suggested at: | |
# http://gis.stackexchange.com/questions/163445/r-solution-for-topologyexception-input-geom-1-is-invalid-self-intersection-er | |
# NE_countries <- gBuffer(NE_countries, byid=TRUE, width=0) | |
# intersecting line with country polygons | |
line.gInt <- gIntersection(split.line, NE_countries) | |
# create a very thin polygon (buffer) out of the intersecting "split line" | |
bf <- gBuffer(line.gInt, byid=TRUE, width=0.000001) | |
# split country polygons using intersecting thin polygon (buffer) | |
NE_countries.split <- gDifference(NE_countries, bf, byid=TRUE) | |
# plot(NE_countries.split) # check map | |
class(NE_countries.split) # is a SpatialPolygons object | |
# ~~~~~~~~~~~ Create graticules ~~~~~~~~~~~ # | |
# create a bounding box - world extent | |
b.box <- as(raster::extent(180, -180, -90, 90), "SpatialPolygons") | |
# assign CRS to box | |
proj4string(b.box) <- WGS84 | |
# create graticules/grid lines from box | |
grid <- gridlines(b.box, | |
easts = seq(from=-180, to=180, by=20), | |
norths = seq(from=-90, to=90, by=10)) | |
# create labels for graticules | |
grid.lbl <- labels(grid, side = 1:4) | |
# transform labels from SpatialPointsDataFrame to a data table that ggplot can use | |
grid.lbl.DT <- data.table(grid.lbl@coords, grid.lbl@data) | |
# prepare labels with regular expression: | |
# - delete unwanted labels | |
grid.lbl.DT[, labels := gsub(pattern="180\\*degree|90\\*degree\\*N|90\\*degree\\*S", replacement="", x=labels)] | |
# - replace pattern "*degree" with "°" (* needs to be escaped with \\) | |
grid.lbl.DT[, lbl := gsub(pattern="\\*degree", replacement="°", x=labels)] | |
# - delete any remaining "*" | |
grid.lbl.DT[, lbl := gsub(pattern="*\\*", replacement="", x=lbl)] | |
# adjust coordinates of labels so that they fit inside the globe | |
grid.lbl.DT[, long := ifelse(coords.x1 %in% c(-180,180), coords.x1*175/180, coords.x1)] | |
grid.lbl.DT[, lat := ifelse(coords.x2 %in% c(-90,90), coords.x2*82/90, coords.x2)] | |
# ~~~~~~~~~~~ Prepare data for ggplot, shift & project coordinates ~~~~~~~~~~~ # | |
# give the PORJ.4 string for Eckert IV projection | |
PROJ <- "+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" | |
# transform graticules from SpatialLines to a data table that ggplot can use | |
grid.DT <- data.table(map_data(SpatialLinesDataFrame(sl=grid, | |
data=data.frame(1:length(grid)), | |
match.ID = FALSE))) | |
# project coordinates | |
# assign matrix of projected coordinates as two columns in data table | |
grid.DT[, c("X","Y") := data.table(project(cbind(long, lat), proj=PROJ))] | |
# project coordinates of labels | |
grid.lbl.DT[, c("X","Y") := data.table(project(cbind(long, lat), proj=PROJ))] | |
# transform split country polygons in a data table that ggplot can use | |
Country.DT <- data.table(map_data(as(NE_countries.split, "SpatialPolygonsDataFrame"))) | |
# Shift coordinates | |
Country.DT[, long.new := long + shift] | |
Country.DT[, long.new := ifelse(long.new > 180, long.new-360, long.new)] | |
# project coordinates | |
Country.DT[, c("X","Y") := data.table(project(cbind(long.new, lat), proj=PROJ))] | |
# ~~~~~~~~~~~ Plot map ~~~~~~~~~~~ # | |
ggplot() + | |
# add projected countries | |
geom_polygon(data = Country.DT, | |
aes(x = X, y = Y, group = group), | |
colour = "gray70", | |
fill = "gray90", | |
size = 0.25) + | |
# add graticules | |
geom_path(data = grid.DT, | |
aes(x = X, y = Y, group = group), | |
linetype = "dotted", colour = "grey50", size = .25) + | |
# add a bounding box (select graticules at edges) | |
geom_path(data = grid.DT[(long %in% c(-180,180) & region == "NS") | |
|(long %in% c(-180,180) & lat %in% c(-90,90) & region == "EW")], | |
aes(x = X, y = Y, group = group), | |
linetype = "solid", colour = "black", size = .3) + | |
# add graticule labels | |
geom_text(data = grid.lbl.DT, # latitude | |
aes(x = X, y = Y, label = lbl), | |
colour = "grey50", size = 2) + | |
# ensures that one unit on the x-axis is the same length as one unit on the y-axis | |
coord_equal() + # same as coord_fixed(ratio = 1) | |
# set empty theme | |
theme_void() | |
# ~~~~~~~~~~~ Save to pdf and png file ~~~~~~~~~~~ # | |
ggsave("World_map_shift.pdf", width=29.7, height=14, units="cm") | |
ggsave("World_map_shift.png", width=29.7, height=14, units="cm", dpi=300) |
Hi , sorry, I'll not have time to look into this. Please check https://www.naturalearthdata.com/downloads/ and see the options or download them manually from there.
Hi there, when using this code I'm having issues with the labels on the longitude lines. The map is Pacific-centred, but the longitudinal line in the centre is now labelled as 0 degrees. For example, New Zealand is shown as being at 20 degrees east, when it should be at 166-170 degrees east. Which part of your code do I need to change to fix the labels? Any help hugely appreciated, thanks!
@corch884 , I do not have time to maintain this anymore (my free time is limited), and it looks like many of the packages will be deprecated soon. Maintaining code is almost a full time job in a ever changing world.
At a first look, I would check the coordinates in grid.lbl.DT
or adjust the label values there.
Hi,
the link to download countries data is not working. Can you provide an updated link?
Thanks