Plotting a certain range in azimuthal projection using ggplot2

Question

I would like to limit a map to a particular range of values (i.e. northern hemisphere) under an azimuthal or orthographic projection in ggplot2. When I set ylim, I receive an error message that relates to geom_path() and certain portions of the map are missing (note the triangle cut out of the United States). Any ideas on how to plot a particular region using a non-cartesian projection without losing portions of the shapefile?

world<-map_data('world') 

p<-ggplot(data=world,aes(x=long,y=lat,group=group)) +
ylim(30,90) +
geom_polygon(fill='skyblue') +
geom_path(color='gray60',aes(x=long,y=lat),size=0.1) +
coord_map("orthographic")
p

Warning message:
Removed 276 rows containing missing values (geom_path). 

Answer 1

I'm not sure if ggplot depends on proj.4, but I solved similar problems in QGIS by defining the projection on a sphere instead an ellipsoid:

Manipulating Azimuthal Equidistant Projections in QGIS

For orthographic projections, it might be better to crop the polygons to the visible hemisphere:

Where did the polygons go after projecting a map in QGIS?

See also: http://docs.ggplot2.org/current/coord_map.html (the last remark).

And here you may find some hints on how to use proj strings: use proj4 to specify Robinson projection with R ggmap and ggplot2 packages?

Here are some more tips:

https://stackoverflow.com/questions/11490493/preventing-partial-mapping-of-coastlines

https://stackoverflow.com/questions/20143851/ggmap-is-unable-to-plot-specific-region

Answer 2

What I suggest below is to generate a polygon with the area of interest and clip your polygons.

library(ggplot2)
library(sp)
library(rgeos)
library(maptools)

data(wrld_simpl) # use this since it is already a SpatialPolygons object

# make a bounding box for the northern hemisphere
g <- gridlines(wrld_simpl, easts=0, norths=30, ndiscr = 200)
x <- c(coordinates(g)[[1]][[1]][,1], rev(coordinates(g)[[1]][[1]][,1]), coordinates(g)[[1]][[1]][1])
y <- c(rep(90, 200), rep(30, 200), 90)
bnds <- cbind(x=x, y=y)
SP <- SpatialPolygons(list(Polygons(list(Polygon(bnds)), "1")), proj4string=CRS(proj4string(wrld_simpl)))

# the code snippet below is courtesy of Roger Bivand https://stat.ethz.ch/pipermail/r-sig-geo/2012-June/015340.html

gI <- gIntersects(wrld_simpl, SP, byid=TRUE)
out <- vector(mode="list", length=length(which(gI)))
ii <- 1
for (i in seq(along=gI)) if (gI[i]) {
  out[[ii]] <- gIntersection(wrld_simpl[i,], SP)
  row.names(out[[ii]]) <- row.names(wrld_simpl)[i]; ii <- ii+1
}
out1 <- do.call("rbind", out)

# fortify SP object to plot wiht ggplot
mymap <- fortify(out1, region='id') 
mymap <- mymap[order(mymap$order), ]

ggplot(data=mymap,aes(x=long,y=lat,group=group)) +
  geom_polygon(fill='skyblue') +
  geom_path(color='gray60',aes(x=long,y=lat),size=0.1) +
  coord_map("orthographic")