Last active
July 12, 2023 08:25
-
-
Save valentinitnelav/c7598fcfc8e53658f66feea9d3bafb40 to your computer and use it in GitHub Desktop.
Pacific centered world map with ggplot (change central/prime meridian)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# ====================================================================================== | |
# 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) |
@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.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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!