Overlap image plot on a Google Map background in R

Question

I'm trying to add this plot of a function defined on Veneto (italian region)

obtained by an image and contour:

image(X,Y,evalmati,col=heat.colors(100), xlab="", ylab="", asp=1,zlim=zlimits,main=title)
contour(X,Y,evalmati,add=T)

(here you can find objects: https://dl.dropboxusercontent.com/u/47720440/bounty.RData)

on a Google Map background.

I tried two ways:

PACKAGE RGoogleMaps

I downloaded the map mbackground

MapVeneto<-GetMap.bbox(lonR=c(10.53,13.18),latR=c(44.7,46.76),size = c(640,640),MINIMUMSIZE=TRUE)
PlotOnStaticMap(MapVeneto)

but i don't know the commands useful to add the plot defined by image and contour to the map

PACKAGE loa

I tried this way:

lat.loa<-NULL
lon.loa<-NULL
z.loa<-NULL
nx=dim(evalmati)[1]
ny=dim(evalmati)[2]
for (i in 1:nx)
{
    for (j in 1:ny)
    {
        if(!is.na(evalmati[i,j]))
        {
            lon.loa<-c(lon.loa,X[i])
            lat.loa<-c(lat.loa,Y[j])
            z.loa<-c(z.loa,evalmati[i,j])
        }
    }
}

GoogleMap(z.loa ~ lat.loa*lon.loa,col.regions=c("red","yellow"),labels=TRUE,contour=TRUE,alpha.regions=list(alpha=.5, alpha=.5),panel=panel.contourplot)

but the plot wasn't like the first one:

in the legend of this plot I have 7 colors, and the plot use only these values. image plot is more accurate.

How can I add image plot to GoogleMaps background?

Answer 1

If the use of a GoogleMap map is not mandatory (e.g. if you only need to visualize the coastline + some depth/altitude information on the map), you could use the package marmap to do what you want. Please note that you will need to install the latest development version of marmap available on github to use readGEBCO.bathy() since the format of the files generated when downloading GEBCO files has been altered recently. The data from the NOAA servers is fine but not very accurate in your region of interest (only one minute resolution vs half a minute for GEBCO). Here is the data from GEBCO I used to produce the map : GEBCO file

library(marmap)

# Get hypsometric and bathymetric data from either NOAA or GEBCO servers
# bath <- getNOAA.bathy(lon1=10, lon2=14, lat1=44, lat2=47, res=1, keep=TRUE)
bath <- readGEBCO.bathy("GEBCO_2014_2D_10.0_44.0_14.0_47.0.nc")

# Create color palettes for sea and land
blues <- c("lightsteelblue4", "lightsteelblue3", "lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))

# Plot the hypsometric/bathymetric map
plot(bath, land=T, im=T, lwd=.03, bpal = list(c(0, max(bath), greys), c(min(bath), 0, blues)))
plot(bath, n=1, add=T, lwd=.5)  # Add coastline

# Transform your data into a bathy object
rownames(evalmati) <- X
colnames(evalmati) <- Y
class(evalmati) <- "bathy"

# Overlay evalmati on the map
plot(evalmati, land=T, im=T, lwd=.1, bpal=col2alpha(heat.colors(100),.7), add=T, drawlabels=TRUE) # use deep= shallow= step= to adjust contour lines
plot(outline.buffer(evalmati),add=TRUE, n=1) # Outline of the data

# Add cities locations and names
library(maps)
map.cities(country="Italy", label=T, minpop=50000)

Since your evalmati data is now a bathy object, you can adjust its appearance on the map like you would for the map background (adjust the number and width of contour lines, adjust the color gradient, etc). plot.bath() uses both image() and contour() so you should be able to get the same results as when you plot with image(). Please take a look at the help for plot.bathy() and the package vignettes for more examples.

Answer 2

I am not realy inside the subject, but Lovelace, R. "Introduction to visualising spatial data in R" might help you https://github.com/Robinlovelace/Creating-maps-in-R/raw/master/intro-spatial-rl.pdf From section "Adding base maps to ggplot2 with ggmap" with small changes and data from https://github.com/Robinlovelace/Creating-maps-in-R/archive/master.zip

library(dplyr)
library(ggmap)
library(rgdal)   

lnd_sport_wgs84 <- readOGR(dsn = "./Creating-maps-in-R-master/data", 
                           layer = "london_sport") %>%
  spTransform(CRS("+init=epsg:4326")) 

lnd_wgs84_f <- lnd_sport_wgs84 %>%
  fortify(region = "ons_label") %>%
  left_join(lnd_sport_wgs84@data, 
            by = c("id" = "ons_label"))

ggmap(get_map(location = bbox(lnd_sport_wgs84) )) + 
  geom_polygon(data = lnd_wgs84_f,
               aes(x = long, y = lat, group = group, fill = Partic_Per),
               alpha = 0.5)