Dear all,

I'm a relatively new with R. However, I realized from colleagues the great
capabilities of R or analyzing geo data.
I'm interested in performing an ordinary kriging to interpolate z value on
a world map sarting from a dataset that (unfortunately) includes rather few
observations (less than 10 location worldwide). Of course, it represents a
somewhat exploratory analysis and also a chance to learn the related
function in R.
my problems are the following:

1) the obtained plot extend the kriging to seas and i would like to "color"
only lands
2)I can't find a way to illustrate the kriging in a more "wave (or bands)
of colors rather that in a extremely smooth way as you see in the attached
figure


My dataset (value.txt) is a .txt file:

Lat,Long,value
60,-95.25,-0.021725
-35.3,131.1499998,-0.020675
48.20833,16.373064,0.0098
19.733333,96.2,0.003978
30.058056,31.228889,0.027548889
64.166667,-51.75,-0.014971667
28.6,77.2,0.00255125
-6.2,106.8,0.01832
46.616667,10.07,-0.0055325
-41.283333,174.766667,0.02026
-20,-56,-0.04445
55.751667,37.617778,-0.026996
-33.9264,18.4227,0.00696

below my commands

  library(maptools);library(fields);library(maps);library(geoR)
#here I import a .shp file representing world coastlines downloaded from
http://www.soest.hawaii.edu/pwessel/gshhs/index.html
  world<-readShapeSpatial("10m_coastline.shp")
  dati<-read.table("value.txt",header=TRUE)
  attach (dati)
  s<-cbind(Long,Lat);
  ml <- likfit(data=dati[,3],coords=s,
               fix.nugget=F,cov.model="exponential",
               ini = c(30, 5),nugget=5);
  sp1<-seq(min(s[,1]),max(s[,1]),length=500);
  sp2<-seq(min(s[,2]),max(s[,2]),length=500);
  sp<-expand.grid(sp1,sp2);
  inCA<-map.where("state",x=sp[,1],y=sp[,2]);
  inCA<-inCA=="world"
  pred<-krige.conv(data=dati[,3],coords=s,locations=sp,
                   krige=krige.control(cov.model="exponential",
                                       cov.pars=c(26.19,13.7),
                                       nugget=4.5763));
  image.plot(sp1,sp2,matrix(pred$predict,500,500),zlim=range(dati[,3]))
  map("world",add=T)

	[[alternative HTML version deleted]]

