Eu tenho dados de levantamento da riqueza de espécies que foram tiradas em vários locais na Baía de Chesapeake, EUA, e gostaria de apresentar graficamente os dados como um "mapa de calor".

Eu tenho um dataframe de lat / long coordenadas e valores de riqueza, que eu convertido em umSpatialPointsDataFrame e usou oautoKrige() função do pacote automap para gerar os valores interpolados.

Primeiro, alguém pode comentar se estou implementando corretamente oautoKrige() função?

Em segundo lugar, estou tendo problemas para plotar os dados e sobrepor um mapa da região. Alternativamente, eu poderia especificar a grade de interpolação para refletir as bordas da Baía (como sugeridoAqui) Alguma idéia de como eu poderia fazer isso e onde poderia obter essa informação? Fornecendo a grade paraautoKrige() parece bastante fácil.

EDIT: Graças a Paul por seu post super útil! Aqui está o que tenho agora. Tendo problemas em fazer com que o ggplot aceite tanto os dados interpolados quanto a projeção do mapa:

<code>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")
</code>