[R-sig-Geo] question about interpreting timelag in output of variogramST

[Rawling, Geoffrey]

Hello all

I have done quite a bit of searching including in the archives of this list
and I cannot find an answer to this question. The script below calculates a
sample spacetime variogram from water levels measured in wells over about
40 years. The variogram plots look fine but I dont understand the values in
the timelag field of the output, it comes out as decimals, with a max value
of 0.52. If I have formatted the input data correctly shouldn't the time
lags be counted in months as the time values of the data points are "month
year", just as the lags for the spatial coordinates are in meters?

Do the decimal values in the timefield of the variogram ST output mean the
time lag as a fraction of the entire time series of data?

I have also formatted the input data as the full date of measurement,
including the day, the time lags then result in days (I think) but with a
maximum of about 190. But the total time series of data spans about 40
years. I would like to use the temporal data spanning a couple of tens of
years to inform the sample variogram. I have tried various values for the
tlags parameter in variogram ST but that has only furthered my confusion.

I am probably missing something obvious about how variogramST works and/or
how to format the temporal data correctly for input to the function.

Any advice is appreciated.

Geoff Rawling

dropbox link to two csv files for script
https://www.dropbox.com/s/w2vcelj4sr5d779/doi_wells.zip?dl=0


# script to calculate experimental spacetime variogram

library(sp)
library(spacetime)
library(gstat)
library(rgdal)

library(ggplot2)
library(xts)

# read in data of interest, doi

#read in data comma separated text files of well locs and wls
doi = read.table("doi.csv", header = TRUE, sep = ",",  stringsAsFactors =
FALSE)

wells = read.table("wells.csv", header = TRUE, sep = ",",  stringsAsFactors
= FALSE)

# basic heat map of wl meas pts thru time

d <- ggplot(doi, aes(mnth_ct, Well_SID )) +
scale_fill_gradient(low="black", high="red")

d + geom_bin2d(binwidth = c(3, 1))

# simplify to create index matrix for input STSDF, each unique well site
gets a number, each unique time of meas gets a number
# n x 2 matrix, each row is a wl meas, 1st column is well site number,
second column is meas time number

index <- as.matrix(doi[,c("Well_SID","T_SID")])

# sort the well locs by increasing  x, then y

wells <- wells[order(wells$LCC_East, wells$LCC_North),]

# plot to check data locations, symbolize w well id
ggplot(doi, aes(LCC_East,LCC_North)) + geom_point()
ggplot(doi, aes(LCC_East,LCC_North)) + geom_text(aes(label=Well_SID))

rownames(wells) <- wells[,1]
wells <- wells[,c(2,3)]


wells = SpatialPoints(wells, proj4string=CRS("+proj=lcc +lat_0=0
+lon_0=-106 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0 +datum=WGS84 +units=m
+no_defs "))

# time of water level measurements
# DateL is from lubridate

time <- as.Date(doi[,"DateL"])

# simplify time to month and year
timeYM <- as.yearmon(time)
str(timeYM)
summary(timeYM)


# PCA trend surface residuals at spatial points and time of measurement

tsr <- doi[,c("PCA.tsR","Point_ID")]

# make spacetime object STSDF

# STSDF(sp, time, data, index, endTime = delta(time))
SO <- STSDF(wells, time, tsr, index)

timeSO <- SO using time
str(timeSO)
str(SO)


# calc exp st variogram

exp_st_vgm <- variogramST(PCA.tsR~1, SO, width = 2000, cutoff = 20000,
assumeRegular=FALSE, na.omit = TRUE)

plot(exp_st_vgm)# plot variograms
plot(exp_st_vgm, wireframe = T) # plot variograms
plot(exp_st_vgm, map=F) # plot variograms

-- 
*********************************************

Geoffrey Rawling

Senior Field Geologist
NM Bureau of Geology and Mineral Resources

Geoffrey.Rawling using nmt.edu

1015 Tijeras Avenue NW, Suite 200
Albuquerque, NM, 87102

505-366-2535 phone
505-999-8307 cell

*********************************************

	[[alternative HTML version deleted]]