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 acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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