Plotting contours on an irregular grid

Question

I have gone through pages and pages of contour plots in R (including many hints on stackoverflow) without success. Here is my data to contour, including adding a map of Rwanda (the data consists of 14 values of longitude, latitude and rain as x,y and z):

Lon Lat Rain 28.92   -2.47   83.4 29.02   -2.68   144 29.25   -1.67   134.7 29.42   -2.07   174.9 29.55   -1.58   151.5 29.57   -2.48   224.1 29.6    -1.5    254.3 29.72   -2.18   173.9 30.03   -1.95   154.8 30.05   -1.6    152.2 30.13   -1.97   126.2 30.33   -1.3    98.5 30.45   -1.81   145.5 30.5    -2.15   151.3 

Here is the code I tried from stackoverflow:

datr <- read.table("Apr0130precip.txt",header=TRUE,sep=",") x <- datr$x y <- datr$y z <- datr$z  require(akima)  fld <- interp(x,y,z)  par(mar=c(5,5,1,1)) filled.contour(fld) 

The interpolation fails.help will be appreciated.

Answer 1

Here are some different possibilites using base R graphics and ggplot. Both simple contours plots, and plots on top of maps are generated.

---

Interpolation

library(akima) fld <- with(df, interp(x = Lon, y = Lat, z = Rain)) 

---

base R plot using filled.contour

filled.contour(x = fld$x,                y = fld$y,                z = fld$z,                color.palette =                  colorRampPalette(c("white", "blue")),                xlab = "Longitude",                ylab = "Latitude",                main = "Rwandan rainfall",                key.title = title(main = "Rain (mm)", cex.main = 1)) 

---

Basic ggplot alternative using geom_tile and stat_contour

library(ggplot2) library(reshape2)  # prepare data in long format df <- melt(fld$z, na.rm = TRUE) names(df) <- c("x", "y", "Rain") df$Lon <- fld$x[df$x] df$Lat <- fld$y[df$y]  ggplot(data = df, aes(x = Lon, y = Lat, z = Rain)) +   geom_tile(aes(fill = Rain)) +   stat_contour() +   ggtitle("Rwandan rainfall") +   xlab("Longitude") +   ylab("Latitude") +   scale_fill_continuous(name = "Rain (mm)",                         low = "white", high = "blue") +   theme(plot.title = element_text(size = 25, face = "bold"),         legend.title = element_text(size = 15),         axis.text = element_text(size = 15),         axis.title.x = element_text(size = 20, vjust = -0.5),         axis.title.y = element_text(size = 20, vjust = 0.2),         legend.text = element_text(size = 10)) 

---

ggplot on a Google Map created by ggmap

# grab a map. get_map creates a raster object library(ggmap) rwanda1 <- get_map(location = c(lon = 29.75, lat = -2),                   zoom = 9,                   maptype = "toner",                   source = "stamen") # alternative map # rwanda2 <- get_map(location = c(lon = 29.75, lat = -2), #                   zoom = 9, #                   maptype = "terrain")  # plot the raster map g1 <- ggmap(rwanda1) g1  # plot map and rain data # use coord_map with default mercator projection g1 +    geom_tile(data = df, aes(x = Lon, y = Lat, z = Rain, fill = Rain), alpha = 0.8) +   stat_contour(data = df, aes(x = Lon, y = Lat, z = Rain)) +   ggtitle("Rwandan rainfall") +   xlab("Longitude") +   ylab("Latitude") +   scale_fill_continuous(name = "Rain (mm)",                         low = "white", high = "blue") +   theme(plot.title = element_text(size = 25, face = "bold"),         legend.title = element_text(size = 15),         axis.text = element_text(size = 15),         axis.title.x = element_text(size = 20, vjust = -0.5),         axis.title.y = element_text(size = 20, vjust = 0.2),         legend.text = element_text(size = 10)) +   coord_map() 

---

ggplot on a map created from shapefile

# Since I don't have your map object, I do like this instead: # get map data from # http://biogeo.ucdavis.edu/data/diva/adm/RWA_adm.zip # unzip files to folder named "rwanda"  # read shapefile with rgdal::readOGR # just try the first out of three shapefiles, which seemed to work. # 'dsn' (data source name) is the folder where the shapefile is located # 'layer' is the name of the shapefile without the .shp extension.  library(rgdal) rwa <- readOGR(dsn = "rwanda", layer = "RWA_adm0") class(rwa) # [1] "SpatialPolygonsDataFrame"  # convert SpatialPolygonsDataFrame object to data.frame rwa2 <- fortify(rwa) class(rwa2) # [1] "data.frame"  # plot map and raindata   ggplot() +    geom_polygon(data = rwa2, aes(x = long, y = lat, group = group),                colour = "black", size = 0.5, fill = "white") +   geom_tile(data = df, aes(x = Lon, y = Lat, z = Rain, fill = Rain), alpha = 0.8) +   stat_contour(data = df, aes(x = Lon, y = Lat, z = Rain)) +   ggtitle("Rwandan rainfall") +   xlab("Longitude") +   ylab("Latitude") +   scale_fill_continuous(name = "Rain (mm)",                         low = "white", high = "blue") +   theme_bw() +   theme(plot.title = element_text(size = 25, face = "bold"),         legend.title = element_text(size = 15),         axis.text = element_text(size = 15),         axis.title.x = element_text(size = 20, vjust = -0.5),         axis.title.y = element_text(size = 20, vjust = 0.2),         legend.text = element_text(size = 10)) +   coord_map() 

---

The interpolation and plotting of your rainfall data could of course be done in a much more sophisticated way, using the nice tools for spatial data in R. Consider my answer a fairly quick and easy start.