[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: creating a variogram with gstat in R
Sorry to be dense about this but I'm still having trouble working out
how gstat operates within R in terms of data import (the data import
phase isn't in the examples in the demos).
OK -- maybe I aimed more at people familiar with R than
people familiar with gstat-standalone, sorry for that.
I import my data fine - and have been trying it with two approaches,
one with the original ascii grid (dataA) and the other with a csv list
of x,y, and value (dataB).
> dataA <- read.table(file="testData.txt", sep=" ", na.strings =
"-9999", skip = 6)
I wouldn't expect success here: asciigrid files only contain z values,
not x or y values (as they're implicitly defined by row,col, grid origin
and cell size values).
If you want to create a grid (without missing values), explore the
S function expand.grid, e.g.
You could get the wanted table by:
predictionlocs = data.frame(expand.grid(x=..., y=...), dataA)
which glues x and y (with ... to be filled in correctly) to the
attribute values (and possible NA's) in dataA.
> dataB <- read.table("testData2.csv", sep=",", header=TRUE, row.names=1)
I hade more success when using read.cvs, but maybe this works.
I don't think the row.names=1 does anything.
In attempting to plot a basic variogram, I do the following and get
the an error message I can't interpret (nor is there much help on line
with this error message):
> v <- variogram(log(A1)~1,loc= ~x+y, dataB)
Error in vector("double", length) : negative length vectors are not
I'm also wondering how the variogram function knows which column is x
and which y (as is specified in the command files), or is the x and y
This is in the manual pages. The loc=~x+y tells the variogram
function that locations are in columns named "x" and "y" in dataB.
With dataA (the asciigrid) I'm not even sure how to create a variogram
as it doesn't have a column header.
Maybe you can export the grid from ArcGIS in a column format. In R,
it is certainly possible, e.g. using the expand.grid() function, but a
more direct way (a function reading the asciigrid header, calling
expand.grid directly with the right parameters) I wouldn't know.
You could post a question on this subject on R-sig-geo, see
Thanks in advance, and apologies for the simplistic questions.
No, thank you for the questions, from which we, developers,
learn from you, users, where you get stuck.