I am trying to produce a high quality vector based map (in greyshades) in order to report locations.

I tried this by getting maps from GADM and SRTM to make a raster as a base for plotting in ggplot2. However, the resulting data frame is too large for plotting. Quesion 1: How to simplify the data frame for ggplot2 without while still obtaining a high quality resolution map?

Here is what I tried for south New Zealand:


nz1 <- getData('GADM', country='NZL', level=1)
nz1 <- subset(nz1,NAME_1 %in% c("Southland","Otago","West Coast"))

nz1c <- gCentroid(nz1) %>% coordinates()

dem1 <- getData("SRTM",lat=nz1c[2],lon=nz1c[1],path=datadir)
dem2 <- getData("SRTM", lat = -45.516667, lon = 168.566667,path=datadir) # City of Athol coordinates
dem3 <- getData("SRTM",lat = -45.866667, lon = 170.5,path=datadir) # City of Dunedin coordinates
dem4 <- getData("SRTM",lat = -44.383333, lon = 168.716667,path=datadir) # Mount Aspiring GPS data -44.383333, 168.716667

dem <- merge(dem1,dem2,dem3,dem4)

dem <- crop(dem,nz1,filename=file.path(datadir,"dem_nz1.tif")

dem.p  <-  rasterToPoints(dem)
df <-  data.frame(dem.p)
colnames(df) = c("lon", "lat", "alt")

p1 <- ggplot(df, aes(lon,lat)) +
  geom_raster(aes(fill = alt))

The data frame has about 25 Million observations Which produces something like this (~110MB PDF):~110MB PDF

This looks already very promising but it is way too heavy.

How could I possibly reduce the data frame or use any other approach to obtain a nice high scalable vector image for my ggplot2 maps?

I just started working with maps in R and basically I just want to produce nice terrain/topographic maps, and I don't want to use ggmaps or similar for this.

  • What is za? nz1? Also datadir, nz0c are missing. I would find the right resolution of the entire grid for you final map. After the crop, try dem1 <- aggregate(dem, fact = 16) and then replot with dem1, then reduce (or increase) "fact" and try the plot again until you find the level that's right. You can use plot(dem, maxpixels = ncell(dem) / (fact * fact)) to get a faster result without having to create points.
    – mdsumner
    Commented Jan 11, 2017 at 4:29
  • Thank you for your comment. Yes, za should be nz1 and nz0c should be nz1c. Sorry for bad example code. I was just trying to show how I got the coordinates for the STRM rasters. At least I corrected that for now but as I am not able to test this and your code now, I report back tomorrow. Thanks.
    – Tony
    Commented Jan 11, 2017 at 10:10
  • 1
    I just can't get the downloads to work, can you list the specific links to those files? getData prints out the links when it dowloads. I don't know what the problem is, getData hasn't worked for me for some time. (Also, your code: datadir is still missing, and rasterVis is not used. )
    – mdsumner
    Commented Jan 11, 2017 at 23:36
  • I already tried your suggestion today and it worked quite well. I found fact = 5 to be a suitable compromise. So thanks anyways, that was a great help. However, didn't had time for reworking the script here. I'll change it, hopefully soon, with file links.
    – Tony
    Commented Jan 12, 2017 at 9:38
  • No worries, great!
    – mdsumner
    Commented Jan 12, 2017 at 11:53

2 Answers 2


Have you tried lowering the resolution of your raster data before creating the plot? Unfortunately, R is still slow for plotting large images so there is a trade-off between resolution and plotting speed.

    # reduce raster resolution
      dem_lower_res <- aggregate(dem, fact=10)

then plot

dem.p  <-  rasterToPoints(dem_lower_res )
df <-  data.frame(dem.p)
colnames(df) = c("lon", "lat", "alt")

p1 <- ggplot(df, aes(lon,lat)) +
  geom_raster(aes(fill = alt))

Here's a way, efficiency is enabled by the new terrainmeshr package, and we convert to gg-ready form with anglr.

Just replace dem with your data, I couldn't get getData to work so I grabbed some AWS tiles.

#> Loading required package: sp
  # whatever you do to get your dem
  el <- cc_location(extent(166, 175, -48, -40), 
                  type = "elevation-tiles-prod", max_tiles = 36)
 dem <- projectRaster(el, crs =  "+proj=longlat") 
 dem[dem<1] <- NA

#> Preparing to download: 20 tiles at zoom = 7 from 
#> https://s3.amazonaws.com/elevation-tiles-prod/geotiff/

mesh <- as.mesh3d(dem)
## plot3d(mesh)
## fortify ggplot()  
fortify.mesh3d <- function(x, ...) {
  idx <- if (!is.null(x$it)) x$it else x$ib
  nc <- dim(idx)[2L]
  idx <- as.vector(idx)
  xx <- x ## workaround the tibble name-steal
  tibble::tibble(x = xx$vb[1L, idx],
                 y = xx$vb[2L, idx],
                 z = xx$vb[3L, idx],
                 group = rep(seq_len(nc), each = 3L))
ggplot(mesh) + geom_polygon(aes(x, y, group = group, fill = z), colour = NA)

Created on 2020-05-15 by the reprex package (v0.3.0)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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