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@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

[hidden email]
1015 Tijeras Avenue NW, Suite 200

Albuquerque, NM, 87102

505-366-2535 phone

505-999-8307 cell

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

[[alternative HTML version deleted]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo