Dear Edzer (et al),
Thanks for the help below, I'm getting closer to pulling this all together.
I thought I would also post the following help about importing an asciigrid
for use in gstat as it proved very useful, which I received from Andy Bunn
on the R-help list . It would be great if this could be made into a gstat
function.
******
If you have exported the grid from Arc using the asciigrid command then
you can read it in with scan or read.table. You can tell R to skip the
six lines of header info and to convert -9999 to NA e.g.,
$ snep.tmin <- read.table(file = "tmin.asc", sep = " ", na.strings =
"-9999", skip = 6)
Check the number of rows and columns to make sure it matches your data
(in Windows, Arc puts a space before the line return at the end of a row
making the resulting R object have one too many columns.)
If that happens, then remove it:
$ snep.tmin <- snep.tmin[,-ncol(snep.tmin)]
(If there is a work around for the read.table command that somebody else
uses then I'd love to hear it.)
For gstat, it looks like it would be helpful to put the grid into a
vector and attach the coordinate information in a data frame? If so,
from the header information take the lower left corner and make your
coordinate columns and join it to the grid data. e.g.,
$ xLLcorner <- -1855500
$ yLLcorner <- -944500
$ cellsize <- 1000
$
$ xURcorner <- xLLcorner + (cellsize * (ncol(snep.tmin) - 1))
$ xLRcorner <- xURcorner
$ xULcorner <- xLLcorner
$
$ yULcorner <- yLLcorner + (cellsize * (nrow(snep.tmin) - 1))
$ yURcorner <- yULcorner
$ yLRcorner <- yLLcorner
$
$ coords <- expand.grid(y = seq(yULcorner, yLRcorner, by = -1000),
+ x = seq(xULcorner, xLRcorner, by = 1000))
$
$ tmin.frame <- data.frame(coords, tmin = as.vector(c(snep.tmin,
recursive = T))) $
$
Watch your signs depends on the coordinate system. From there you can
krige or whatever easily.
----- Original Message -----
From: "Edzer J. Pebesma" <[EMAIL PROTECTED]>
To: "femke" <[EMAIL PROTECTED]>; <[EMAIL PROTECTED]>
Sent: Wednesday, February 18, 2004 3:21 AM
Subject: Re: creating a variogram with gstat in R
> femke wrote:
>
> >
> > Hello again,
> >
> > 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.
>
> expand.grid(x=1:10,y=1:10)
>
> 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
> > allowed
> > 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
> > coded in?
>
> 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
>
> https://www.stat.math.ethz.ch/mailman/listinfo/r-sig-geo
>
> > Please help.
> >
> > 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.
> --
> Edzer
>
