2

I have a raster that I'd like to plot in the Azimuthal Equidistant projection with a grid and labels. So far I've been able to do that. However, I cannot get the labels inside the plotting area. They are always cut off.

I have tried most of the par() options such as playing with par(mar), par(oma), and par(usr). Hard to tell, but I think the par("usr") option does not work with raster::plot. I also tried playing around with the extent of the raster, since R bases the plot size on the coordinate range of the initial data, but that also did not work. Below is my code and cut off labels.

GSmean<-stack("[file][1]") 
parallels <- c(80, 70, 60)
meridians <- seq(-180, 160, by = 30)
mlim <- c(-180, 180)
plim <- c(50, 80)
grat <- graticule(lons = meridians, lats = parallels, xlim = mlim, ylim = plim, proj = projection(GSmean))
labs <- graticule_labels(meridians, parallels, xline = 360, yline =46, proj = projection(GSmean))

#new plot
breaks<-c(-1,0,1,2,999)
cols <- c("darkgreen", "darkred","gray81","gray81")
plot(GSmean, col=cols,breaks=breaks,axes=F,legend=FALSE,box=FALSE)
plot(grat, add = TRUE, lty = 3,col="honeydew3")
text(labs, lab = parse(text= labs$lab), 
     col= c("gray30", "gray30")[labs$islon + 1], cex = 0.85)

enter image description here

2

Try adding xpd = NA to the text() call to prevent clipping.

Here's a stab at an example:

library(raster)
library(graticule)
p <- projectExtent(raster(extent(-180, 180, 50, 90), crs = "+proj=longlat"), 
               "+proj=laea +lat_0=90")
res(p) <- 50000
data("wrld_simpl", package = "maptools")

p <- rasterize(sp::spTransform(subset(wrld_simpl, coordinates(wrld_simpl)[,2] > 0), projection(p)), p)
parallels <- c(80, 70, 60)
meridians <- seq(-180, 160, by = 30)    
mlim <- c(-180, 180)
plim <- c(50, 80)
grat <- graticule(lons = meridians, lats = parallels, xlim = mlim, ylim = plim, proj = projection(p))
labs <- graticule_labels(meridians, parallels, xline = 360, yline =46, proj = projection(p))

#new plot
plot(p, axes=F,legend=FALSE,box=FALSE)
plot(grat, add = TRUE, lty = 3,col="honeydew3")
text(labs, lab = parse(text= labs$lab), 
 col= c("gray30", "gray30")[labs$islon + 1], cex = 0.85,xpd = NA)

enter image description here

  • That did it. I knew it was somewhere little.... Also, @mdsumner, I used your blog post for help with adding the lat/lon lines and labels. Thank you for the work you do. – webbe Jan 8 at 12:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.