Smoothing geom_tile map in ggplot2 - Interpolating data

Question

I'm trying to plot an spatial map from a lon-lat data set. Following the example in https://stackoverflow.com/a/19339663/709777 I am able to almost get the desired result

This is the code used (most of it from the previously cited question)

library(rgdal)
library(ggplot)
# read province limits
CV = readOGR(dsn=".", layer="poligonos_provincia_etrs89")
CV.pr=subset(CV, CV$CODINE == "46" | CV$CODINE == "12" | CV$CODINE == "03" )
# convert object to data.frame
CV2 <- fortify(CV.pr)

# read spatial data
datos.uvi=read.csv("salida-mapa-036.dat",sep=",",header=T,na.strings="-99.9")
uvi.temp<-datos.uvi[,c("longitud","latitud","RGlobal")]
colnames(uvi.temp)<-c("long","lat","RGlobal")

# plot map
ggplot() +
  geom_tile(data = uvi.temp, aes(x = long, y = lat, z = RGlobal, fill = RGlobal), alpha = 0.4) +
  stat_contour(data = uvi.temp, aes(x = long, y = lat, z = RGlobal)) +
  geom_path(data = CV2, aes(x = long, y = lat, group = group),color="black", size=0.6) +
  ggtitle("Previsión UVI - CV") +  xlab("Longitud") +  ylab("Latitud") +
  scale_fill_continuous(name = "UVI", low = "white", high = "red") +
  theme_bw() +
  coord_map()

The point is that I don't want a tile map but a smooth gradient map. This warning appears when running the code:

Warning message: Computation failed in stat_contour(): (list) object cannot be coerced to type 'double'

I assume this is because of the data type in the data frame (three columns with lon,lat and data) shown below

dput(uvi.temp)
structure(list(long = c(-1.25, -1, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1.25, -1, -0.75, -0.5, -0.25, 0, 
0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 
0.75, -1.25, -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, 
-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, 
-0.5, -0.25, 0, 0.25, 0.5, 0.75, -1.25, -1, -0.75, -0.5, -0.25, 
0, 0.25, 0.5, 0.75), lat = c(37.75, 37.75, 37.75, 37.75, 37.75, 
37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 38, 37.75, 37.75, 37.75, 
37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 38, 38, 38, 38, 38, 
38, 38, 38, 38, 38.25, 38.25, 38.25, 38.25, 38.25, 38.25, 38.25, 
38.25, 38.25, 38.5, 38.5, 38.5, 38.5, 38.5, 38.5, 38.5, 38.5, 
38.5, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 
38.75, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39.25, 39.25, 39.25, 
39.25, 39.25, 39.25, 39.25, 39.25, 39.25, 39.5, 39.5, 39.5, 39.5, 
39.5, 39.5, 39.5, 39.5, 39.5, 39.75, 39.75, 39.75, 39.75, 39.75, 
39.75, 39.75, 39.75, 39.75, 40, 40, 40, 40, 40, 40, 40, 40, 40, 
40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 
40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.75, 
40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 37.75, 
37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 37.75, 38, 38, 
38, 38, 38, 38, 38, 38, 38, 38.25, 38.25, 38.25, 38.25, 38.25, 
38.25, 38.25, 38.25, 38.25, 38.5, 38.5, 38.5, 38.5, 38.5, 38.5, 
38.5, 38.5, 38.5, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 38.75, 
38.75, 38.75, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39.25, 39.25, 
39.25, 39.25, 39.25, 39.25, 39.25, 39.25, 39.25, 39.5, 39.5, 
39.5, 39.5, 39.5, 39.5, 39.5, 39.5, 39.5, 39.75, 39.75, 39.75, 
39.75, 39.75, 39.75, 39.75, 39.75, 39.75, 40, 40, 40, 40, 40, 
40, 40, 40, 40, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 40.25, 
40.25, 40.25, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 40.5, 
40.5, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 40.75, 
40.75), RGlobal = c(469.4264, 467.3801, 469.4264, 467.3801, 522.7478, 
443.958, 441.283, 509.0632, 435.3009, 432.4717, 566.5648, 527.1168, 
469.4264, 467.3801, 522.7478, 443.958, 441.283, 509.0632, 435.3009, 
432.4717, 566.5648, 527.1168, 459.3489, 457.781, 439.8998, 437.827, 
435.4037, 433.0625, 430.2155, 564.3985, 592.4016, 455.9498, 455.622, 
443.1096, 431.0866, 500.2498, 428.5327, 427.8865, 561.918, 626.1251, 
456.2038, 458.2242, 459.7855, 431.4439, 492.9174, 425.223, 424.3311, 
490.1372, 450.9101, 449.7869, 586.2458, 455.2244, 452.9787, 514.1946, 
425.5438, 419.4045, 484.3133, 441.4613, 445.6218, 575.9768, 509.4317, 
437.0949, 411.3335, 412.9341, 412.7781, 411.7847, 432.1784, 435.553, 
566.9114, 563.8876, 414.4346, 404.6672, 405.9609, 405.5721, 470.3917, 
426.5352, 427.227, 552.512, 584.1269, 399.0303, 398.1212, 398.3982, 
398.6394, 462.0748, 413.0668, 417.8466, 481.0621, 420.5388, 407.8954, 
518.209, 390.5623, 390.4963, 453.643, 405.8642, 405.2392, 464.2149, 
413.4183, 409.0698, 525.5512, 447.1689, 383.9111, 384.0972, 400.2987, 
403.8821, 405.2475, 405.1649, 402.3291, 519.2192, 517.2463, 379.0551, 
379.1682, 393.9954, 396.8607, 454.2946, 394.8534, 386.8805, 498.197, 
538.8216, 387.8813, 375.4716, 374.1738, 377.1193, 437.4277, 385.8378, 
381.4342, 422.0693, 365.6135, 373.4164, 497.8519, 469.4264, 467.3801, 
522.7478, 443.958, 441.283, 509.0632, 435.3009, 432.4717, 566.5648, 
527.1168, 459.3489, 457.781, 439.8998, 437.827, 435.4037, 433.0625, 
430.2155, 564.3985, 592.4016, 455.9498, 455.622, 443.1096, 431.0866, 
500.2498, 428.5327, 427.8865, 561.918, 626.1251, 456.2038, 458.2242, 
459.7855, 431.4439, 492.9174, 425.223, 424.3311, 490.1372, 450.9101, 
449.7869, 586.2458, 455.2244, 452.9787, 514.1946, 425.5438, 419.4045, 
484.3133, 441.4613, 445.6218, 575.9768, 509.4317, 437.0949, 411.3335, 
412.9341, 412.7781, 411.7847, 432.1784, 435.553, 566.9114, 563.8876, 
414.4346, 404.6672, 405.9609, 405.5721, 470.3917, 426.5352, 427.227, 
552.512, 584.1269, 399.0303, 398.1212, 398.3982, 398.6394, 462.0748, 
413.0668, 417.8466, 481.0621, 420.5388, 407.8954, 518.209, 390.5623, 
390.4963, 453.643, 405.8642, 405.2392, 464.2149, 413.4183, 409.0698, 
525.5512, 447.1689, 383.9111, 384.0972, 400.2987, 403.8821, 405.2475, 
405.1649, 402.3291, 519.2192, 517.2463, 379.0551, 379.1682, 393.9954, 
396.8607, 454.2946, 394.8534, 386.8805, 498.197, 538.8216, 387.8813, 
375.4716, 374.1738, 377.1193, 437.4277, 385.8378, 381.4342, 422.0693, 
365.6135, 373.4164, 497.8519)), .Names = c("long", "lat", "RGlobal"
), class = "data.frame", row.names = c(NA, -246L))

Any idea?

Thanks in advance

Some people are reluctant to download data. Perhaps you can give some code for a reproducible example (perhaps using `dput`). — Axeman, Feb 01 '16 at 14:51
@Axeman You are right, I put the link to shorten the question but maybe this way it will reach more people — pacomet, Feb 02 '16 at 08:07
Do you just want to interpolate your coarse grid data onto a much finer grid so it looks smoother? — Spacedman, Feb 02 '16 at 08:27
You have duplicated data in the _uvi.temp_ data frame. If you remove duplicates (`uvi.temp <- uvi.temp[!duplicated(uvi.temp), ]`), then `stat_contour()` will work. — oshun, Feb 02 '16 at 18:43
Hi @oshun Now contours appear, thanks. I have to check the script that generates data to avoid this. But what I actually need is to smooth the tiles to generate a smooth gradient map. — pacomet, Feb 03 '16 at 08:59

score 3 · Accepted Answer · edited Apr 13 '17 at 12:33

Finally I could find a solution, not a perfect one but a good enough one. I'll keep searching for better interpolation methods. From the question in https://gis.stackexchange.com/q/169184/9227 I made this code

library(ggplot2)
library(gstat)
library(sp)
library(maptools)
library(rgdal)

# Reading three data frames
datos.uvi.1 <- read.csv(file = "./salida-mapa-012.dat",header = TRUE)
datos.uvi.2 <- read.csv(file = "./salida-mapa-036.dat",header = TRUE)
datos.uvi.3 <- read.csv(file = "./salida-mapa-060.dat",header = TRUE)

# Looking for RGlobal max
new.uvi <- data.frame(datos.uvi.1$RGlobal,datos.uvi.2$RGlobal,datos.uvi.3$RGlobal)
max=apply(new.uvi, 1, max, na.rm=FALSE)
datos.uvi=data.frame(datos.uvi.1$longitud,datos.uvi.1$latitud,max)
colnames(datos.uvi)<-c("longitud","latitud","RGlobal")

# Define x & y as longitude and latitude
datos.uvi$x <- datos.uvi$longitud 
datos.uvi$y <- datos.uvi$latitud

coordinates(datos.uvi) = ~x + y

# Reading shapefiles for province limits (black and blue limits in the map)
provincias = readOGR(dsn=".", layer="poligonos_provincia_etrs89")
pr1=subset(provincias, provincias$CODINE == "46" | provincias$CODINE == "12" | provincias$CODINE == "03" )
pr2=subset(provincias, provincias$CODINE == "43" | provincias$CODINE == "16" | provincias$CODINE == "30"  | provincias$CODINE == "02" | provincias$CODINE == "50" )
pr1 <- fortify(pr1)
pr2 <- fortify(pr2)

# Interpolation area
x.range <- as.numeric(c(-1.75, 1))  # min/max longitude
y.range <- as.numeric(c(37.5, 41))  # min/max latitude

# Create a gridded structure
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1], to = y.range[2], by = 0.1))
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE

#Interpolate surface and fix the output. Apply idw model for the data
idw <- idw(formula = RGlobal ~ 1, locations = datos.uvi, newdata = grd)  
idw.output = as.data.frame(idw) 
names(idw.output)[1:3] <- c("long", "lat", "RGlobalmax")

# Plot
  ggplot() + geom_tile(data = idw.output, alpha = 0.8, aes(x = long, y = lat, fill = RGlobalmax)) +
    scale_fill_gradient(low = "cyan", high = "orange",name = "UVI") +
    geom_path(data = pr2, aes(long, lat, group = group), colour = "blue") +
    geom_path(data = pr1, aes(long, lat, group = group), colour = "black") +
    coord_map(xlim = c(-1.7, 1),ylim = c(37.6,40.9)) +
    ggtitle("Previsión UVI - DD/MM/YYYY") + xlab(" ") + ylab(" ")

And the final output

Smoothing geom_tile map in ggplot2 - Interpolating data

1 Answers1

Linked