Нанесение интерполированных данных на карту

У меня есть данные обзора видового богатства, которые были получены в различных местах в Чесапикском заливе, США, и я хотел бы графически представить данные в виде «тепловой карты».

У меня есть кадр данных с координатами широты/долготы и значениями насыщенности, которые я преобразовал в SpatialPointsDataFrameи использовал функцию autoKrige()из пакета automap для генерации интерполированных значений.

Во-первых, кто-нибудь может прокомментировать, правильно ли я реализую функцию autoKrige()?

Во-вторых, у меня возникли проблемы с нанесением данных и наложением карты региона. В качестве альтернативы, могу ли я указать сетку интерполяции, чтобы отразить границы залива (, как предлагается здесь)? Любые мысли о том, как я могу это сделать и где я могу получить эту информацию? Задать сетку в autoKrige()достаточно просто.

РЕДАКТИРОВАТЬ:Спасибо Полу за его очень полезный пост! Вот что у меня есть сейчас. Возникли проблемы с тем, чтобы заставить ggplot принимать как интерполированные данные, так и картографическую проекцию:

require(rgdal)
require(automap)
#Generate lat/long coordinates and richness data
set.seed(6)
df=data.frame(
  lat=sample(seq(36.9,39.3,by=0.01),100,rep=T),
  long=sample(seq(-76.5,-76,by=0.01),100,rep=T),
  fd=runif(10,0,10))
initial.df=df

#Convert dataframe into SpatialPointsDataFrame
coordinates(df)=~long+lat

#Project latlong coordinates onto an ellipse
proj4string(df)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
#+proj = the type of projection (lat/long)
#+ellps and +datum = the irregularity in the ellipse represented by planet earth

#Transform the projection into Euclidean distances
project_df=spTransform(df, CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84")) #projInfo(type="proj")

#Perform the interpolation using kriging
kr=autoKrige(fd~1,project_df)
#Extract the output and convert to dataframe for easy plotting with ggplot2
kr.output=as.data.frame(kr$krige_output)
#Plot the output
#Load the map data for the Chesapeake Bay
cb=data.frame(map("state",xlim=range(initial.df$long),ylim=range(initial.df$lat),plot=F)[c("x","y")])

ggplot()+
  geom_tile(data=kr.output,aes(x=x1,y=x2,fill=var1.pred))+  
  geom_path(data=cb,aes(x=x,y=y))+
  coord_map(projection="mercator")