Trouble with universal kriging

Dominik, if you want to receive help, please explain why you are
frustrated, and instead of claims of "won't work", explain in detail
what you did and what happened, preferably with a well-prepared example
that people on the list can reproduce without downloading data first.

Hi, I am trying to follow the examples in  Bivand et al 'Applied Spatial
Data Analysis' and the examples at
http://casoilresource.lawr.ucdavis.edu/drupal/node/442 but I am frustrated.
I finally managed an ordinary kriging example to work but with little
understanding of why/when it should work.
I posted some data: https://www.dropbox.com/s/rtg7gh67ziai3ec/mar2011.RData
Another item I've noticed, if I  define the projection with
projection(object)='+proj=lonlat +datum=NAD83'  for any of my spatial
objects it won't work.  It only works with an defined projection. any help
would be appreciated. Thanks.

load('mar2011.RData')

vt=variogram(snotel~recon,data=mar2011)

vt.fit=fit.variogram(vt,vgm(psill=0.025,nugget=0.015,model='Sph',range=2))

#plot(vt,vt.fit)

ok<-krige(snotel~1,locations=mar2011,newdata=swe.grid,model=vt.fit)

uk<-krige(snotel~mar2011$recon,locations=mar2011,newdata=swe.grid,model=
vt.fit)

#spplot(ok)
Dominik

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Edzer Pebesma
Institute for Geoinformatics (ifgi), University of M?nster
Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251
83 33081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de

Dominik, if you want to receive help, please explain why you are
frustrated, and instead of claims of "won't work", explain in detail
what you did and what happened, preferably with a well-prepared example
that people on the list can reproduce without downloading data first.


On 06/20/2013 09:57 PM, Dominik Schneider wrote:

Hi, I am trying to follow the examples in  Bivand et al 'Applied Spatial
Data Analysis' and the examples at
http://casoilresource.lawr.ucdavis.edu/drupal/node/442 but I am frustrated.
I finally managed an ordinary kriging example to work but with little
understanding of why/when it should work.

Well, if you read chapter 8, you should gain sufficient understanding, 
right? The online examples are not intended to be self-explanatory; they 
all work, because they are tested nightly on a fully updated system; the 
online versions have been modified to keep up with successive changes in 
R and the packages used.

I posted some data: https://www.dropbox.com/s/rtg7gh67ziai3ec/mar2011.RData
Another item I've noticed, if I  define the projection with
projection(object)='+proj=lonlat +datum=NAD83'  for any of my spatial
objects it won't work.  It only works with an defined projection. any help
would be appreciated. Thanks.

Do also pay attention to details:

"+proj=lonlat +datum=NAD83"

is valid PROJ.4, but not recognised in the CRS class, it should be:

"+proj=longlat +datum=NAD83"

All of "latlon, "latlong", "lonlat", and "longlat" are valid PROJ.4, but 
the first two imply northings before eastings, and we decided ten years 
ago to use only "longlat". The documentation is in ?proj4string in sp, but 
the detail is not repeated in ?projection in raster, which does link to 
?"CRS-class" in sp and "spTransform-methods" in rgdal. We will try to make 
sure that use of anything other than "longlat" is an error.

Always include your sessionInfo(), which gives details of your platform 
and the versions of installed packages.

Hope this clarifies,

Roger

load('mar2011.RData')

vt=variogram(snotel~recon,data=mar2011)

vt.fit=fit.variogram(vt,vgm(psill=0.025,nugget=0.015,model='Sph',range=2))

#plot(vt,vt.fit)

ok<-krige(snotel~1,locations=mar2011,newdata=swe.grid,model=vt.fit)

uk<-krige(snotel~mar2011$recon,locations=mar2011,newdata=swe.grid,model=
vt.fit)

#spplot(ok)
Dominik

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no

Thanks for the heads up on +proj. I had not noticed that difference.

I read chapter 8, twice now.  It is not clear to me what the error message
I'm getting means:
Error in gstat.formula.predict(d$formula, newdata, na.action = na.action,  : 
  NROW(locs) != NROW(X): this should not occur
In addition: Warning messages:
1: 'newdata' had 5742 rows but variable(s) found have 200 rows 
2: 'newdata' had 5742 rows but variable(s) found have 200 rows

I thought newdata was supposed to be the new grid on which to predict.

Here is code for working data.
snotel=runif(200,0,3)
recon=sample(snotel,200)
recon[1:50]=snotel[1:50] # i'm not sure how to introduce spatial correlation
in this random data.
x=runif(200,-112.25,-104.125)
y=runif(200,33,43.75)

m2=data.frame(x,y,snotel,recon)
coordinates(m2)=~x+y
projection(m2)='+proj=longlat +datum=NAD83'

swe.grid=expand.grid(x=seq(-112.25,-104.125,0.125),y=seq(33,43.75,0.125))
coordinates(swe.grid) <- ~ x+y
projection(swe.grid)='+proj=longlat +datum=NAD83'
gridded(swe.grid) <- TRUE


#SessionInfo()

sessionInfo()

R version 2.15.2 (2012-10-26)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
base     

other attached packages:
[1] raster_2.1-25   maptools_0.8-23 lattice_0.20-13 foreign_0.8-52 
gstat_1.0-16   
[6] sp_1.0-8       

loaded via a namespace (and not attached):
[1] intervals_0.14.0 spacetime_1.0-4  tools_2.15.2     xts_0.9-3       
[5] zoo_1.7-9   




--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Trouble-with-universal-kriging-tp7583848p7583875.html
Sent from the R-sig-geo mailing list archive at Nabble.com.