[R-sig-Geo] sample semivariogram use, retrieving anisotropic parameters, subsampling

[Paul Hiemstra]

Hi Annette,

If you are interested in automating variogram fitting, take a look at
the autofitVariogram function from the automap package. A small example:

library(automap)
data(meuse)
coordinates(meuse) =~ x+y
variogram = autofitVariogram(zinc~1,meuse)
plot(variogram)

variogram = autofitVariogram(zinc ~ soil + ffreq + dist, meuse)
plot(variogram)

The result of autofitVariogram is both the sample variogram and the
fitted nugget. Extracting elements such as the nugget is quite simple:

nug = variogram$var_model$psill[1]
sill = variogram$var_model$psill[2] + nug
range = variogram$var_model$range[2]

hope this helps,
Paul

On 03/16/2011 03:27 PM, Annette H Elmore wrote:
> Hi,
>
> I'm new to this forum and would like to retrieve anisotropic parameters 
> from a sample semi-variogram.  I don't need to krig, because I'm already 
> working at 1 m resolution and don't have to refine any further.  Instead, 
> I'm just looking for an anisotropic description of the autocorrelation 
> patterns in the area where I'm working.  I have three  questions that I 
> hope you can address, for I've been going in circles with them, a bit:
>
> I think that I don't want to fit a spherical model to the the experimental 
> variogram, because I want to know what's actually coming out of the 
> observations themselves, but I am unclear about how the initial 
> experimental parameters are generated by gstat, prior to using them to 
> fit/converge with a model of choice (for smoothing/kriging).  Perhaps the 
> initial raw "curve" parameters aren't a good representation of the 
> experimental cloud and curve fitting is necessary to best represent my 
> data?  I will be comparing scenes through time, so my most important 
> criteria for generating the nugget, sill, and range, is consistency.  I'm 
> less concerned about smoothing in a way that would permit kriging, 
> particularly if smoothing moves me away from the raw data in a way that 
> might be inconsistent between years.  So if you have any insight into 
> this, I would love to know more.
>
> Is it possible to retrieve a nugget, sill, and range from a sample 
> variogram?  Ideally, I want to take the parameters for each image and add 
> them as the last row in an existing table.  I notice that in the gstat 
> program (non-R version) the user can tweak with an interactive interface, 
> but what I'm really after is a way to automate all of this as much as 
> possible.
>
> Is it possible to include distance tolerances for the point pairs so that 
> you can subsample the data to be evaluated, playing around with using 
> points at varying distance lengths away from one another in the 
> calculations?  For instance, if I have regularly spaced points that are 
> about 1-1.5 m apart, but want to select pairs as if only points 4 - 5 m 
> apart were present, is there a way to do that?  I could do this by 
> fiddling with the input shapefile (deleting intervening points that I 
> didn't want evaluated) before I import it, but it would be much nicer to 
> do it here if it's possible.
>
> If it's of any help, this is the code I have produced, so far:
>  
> library(maptools)
> library(gstat)
> baseCRS<-CRS("+proj=utm+zone=17+datum=NAD83")
> test <- readShapePoints("C:\\89802_test.shp", proj4string = baseCRS, 
> verbose = TRUE)
> summary(test)
> test.v <- variogram(GRID_CODE~coords.x1+coords.x2, data=test,  alpha = 
> c(0,45,90,135), 
>         cutoff =41, width = 2)
>
> Thanks,
> Annie
>  
> * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> Annette H. Elmore
> Land Cover Dynamics & Environmental Processes
> Eastern Geographic Science Center
> U.S. Geological Survey
> 12201 Sunrise Valley Drive
> Reston, VA  20192
> Email:  aelmore at usgs.gov
> Voice:  703 648 4805
> Fax:  703 648 4603
> * * * * * * * * * * * * * *
>
> If I don't know I don't know
>                                                I think I know
>
> If I don't know I know
>                                       I think I don't know
>
> -- Knots
>     R.D. Laing
>
> 	[[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


-- 
Paul Hiemstra, MSc
Global Climate Division
Royal Netherlands Meteorological Institute (KNMI)
Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39
P.O. Box 201 | 3730 AE | De Bilt
tel: +31 30 2206 494

http://intamap.geo.uu.nl/~paul
http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770