0

I am trying to plot a heatmap of a country with some points that are probabilities of occurrence of a event. What I did up to now is next:

library(raster)
library(ggplot2)
Uruguay  <- getData("GADM",country="Uruguay",level=0)
ggplot(Uruguay,aes(x=long,y=lat,group=group)) + 
  ggplot2::lims(x = c(-60, -50), y = c(-35, -30))+
  geom_polygon(aes(x = long, y = lat, group = group, fill=id),color="grey30")+
  coord_map(xlim=c(-1,1)+bbox(Uruguay)["x",],ylim=c(-1,1)+bbox(Uruguay)["y",])+
  scale_fill_discrete(guide="none")+
  theme_bw()+theme(panel.grid=element_blank())

country border

my data to produce the heatmap is

prob <- c(10,20,90,40)
lat <- c(-30.52,-32.04,-33.16,-34.28)
long <- c(-57.40,-55.45,-56.35,-56.40)
data <- data.frame(prob, lat, long)

I think that using ggplot2::stat_density2d and ggplot2::scale_fill_gradientn is the way to go but I don't know how to implement it. I want to produce a heatmap like that

enter image description here

Any help is Welcome. Thanks in advance.

Ronak Shah
  • 377,200
  • 20
  • 156
  • 213

2 Answers2

0

To plot the example data you could just use plot

library(raster)
Uruguay  <- getData("GADM",country="Uruguay",level=0)
plot(Uruguay, col="orange")

As for the map you want to make, there are a lot of choices involved. But here is a basic example

prob <- c(10,20,90,40)
lat <- c(-30.52,-32.04,-33.16,-34.28)
long <- c(-57.40,-55.45,-56.35,-56.40)
data <- data.frame(prob, lat, long)

r <- raster(Uruguay, res=.5)
x <- rasterize(cbind(long, lat), r, prob)
plot(x)
lines(Uruguay)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • but how can span I the colors to the rest of the territory? – Maximiliano Verocai Mar 23 '21 at 20:00
  • You need to explicitly model the surface from the data you have. Kriging, 2D convolution, etc. Then you can plot. – Dan Slone Mar 23 '21 at 20:23
  • I assumed you had many more observations. If you only have three points, what assumptions do you use to fill the gaps? See ?interpolate for technical possibilities, but whether these make sense is a very different type of question (methodological, not coding) – Robert Hijmans Mar 23 '21 at 23:31
0

Finally I could get what I wanted. Henrik's answer in this post was very helpful I share the code with you

library(raster)
library(reshape2)
library(ggplot2)
Uruguay  <- getData("GADM",country="Uruguay",level=1)

#invented data
    prob <- c(5, 90,10,15,99,40,90,25,70,90)
    lat <- c(-31,-31.2,-31.3,-34,-32.5,-32.6,-33.7,-34.9,-34.2,-32.5)
    long <- c(-58.3,-55.1,-57.3,-58.4,-56.5,-54,-57.7,-55.8,-54.1,-53.5)
    prueba <- data.frame(prob, lat, long)
    


library(akima)
fld <- with(prueba, interp(x = long, y = lat, z = prob))

class(Uruguay)
uru <- fortify(Uruguay)

library(reshape2)

# prepare data in long format
df <- melt(fld$z, na.rm = TRUE)
names(df) <- c("x", "y", "prob")
df$long <- fld$x[df$x]
df$lat <- fld$y[df$y]


ggplot() + 
  geom_polygon(data = uru, aes(x = long, y = lat, group = group),
               colour = "black", size = 0.5, fill = "white") +
  geom_tile(data = df, aes(x = long, y = lat, z = prob, fill = prob), alpha = 0.8) +
  ggtitle("Frost probability") +
  xlab("Longitude") +
  ylab("Latitude") +
  scale_fill_continuous(name = "Probability (%)",
                        low = "red", high = "blue") +
  theme_bw() +`enter code here`
  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()