r-spatial / gstat

Spatial and spatio-temporal geostatistical modelling, prediction and simulation
http://r-spatial.github.io/gstat/
GNU General Public License v2.0
194 stars 49 forks source link

variogramST with weekly averages #35

Open andresleslielira opened 5 years ago

andresleslielira commented 5 years ago

Hello Edzer, thank you very much for gstat! I apologize beforehand for my poor english and for my short R experience (I learned it one week ago). I think my question is relatively simple. I want to do kriging using PM10 daily data for 8 static stations in Santiago, Chile, from 1997-2012. I explain what I've done so far, and then i explain my question.

DATA: I have first column with days from 1997 to 2012 and 8 columns with PM10 station data.

1) Using STATA I built weekly averages of PM10 for each station (one column each) 2) When (mean) collapsing by week, I kept the median day of each week so that each row has a date. 3) I imported to R this data, and coded :time <-as.POSIXlt(data$date) with no error

[1] 1997-04-05 UTC [2] 1997-04-12 UTC.... [819] 2012-12-27 UTC

4) I imported stations with it's coordinates and built its projections.

coordinates(stations)=~longitud+latitud proj4string(stations) <- CRS("+proj=longlat + ellps=WGS84")

5) Built the vector of data with 8 stations, ordered: PM10<-as.vector(cbind(data $PM10bosque,data $PM10cerrillos,data $PM10cerronavia,data $PM10condes,data $PM10florida,data $PM10independencia,data $PM10parqueoh,data $PM10pudahuel)) PM10<-data.frame(PM10)

6) Coded: DFPM=STFDF(stations, time, PM10) 7) Then calculated sample variogram this way: varPM10 <- variogramST(PM10~1,data=DFPM,tunit="days",assumeRegular=F,na.omit=T)

Output plot looks like: image

Which looks rather strange actually, so i wondered: QUESTION: Does variogramST (with tunit="days") consider each row of data as one consecutive day? is one able to use it while having one observation per week (a weekly average in this example)?

BenGraeler commented 5 years ago

No, gstat respects the temporal gaps between your data. This can also clearly be seen in your plot, where the "rows" are all 7 days apart (0,7,14,21,... days difference). The argument tunit has only an effect when the ST-structure is a STIDF and controls the unit of the tlags argument. Your plot could be influenced by a recurrent pattern of 3 and 10 weeks

andresleslielira commented 5 years ago

Thank you very much.

andresleslielira commented 5 years ago

Hello again, I am REALLY grateful for your help. Feels great to have support from you guys, considering I am new to this

I am posting a new question here because the context is the same and has to do with how KrigeST does it's prediction. I am testing the whole procedure now with daily data for 2008. I tried two things which gave raise to different questions (in bold).

*FIRST EXPERIMENT

I did the same as before (up to point 6.) but after coding DFPM=STFDF(stations, time, PM10), I coded DFPM<-as(DFPM,"STSDF") because I am working with missing data

Then estimated the variogram and its modelling with: varPM10 <- variogramST(PM10~1,data=DFPM,tunit="days",assumeRegular=F,na.omit=T) sepVgm <- vgmST("separable",space=vgm(0,"Exp", 8, 700),time =vgm(200,"Exp", 15, 700), sill=100) sepVgm <- fit.StVariogram(varPM10, sepVgm)

Which results in: image

Then I used KrigeST this way: gridPM10 <-STF(centroids,time) (centroids defined previously the same way as stations) krigedPM10<-krigeST(PM10~1, DFPM, newdata=gridPM10,modelList=sepVgm)

The result of ploting one station data and kriged data for that station county's centroid is: image

, which seems as if the estimation ocurrs for time windows by a set of dates. Does this stairwise strange look have to do with something of the function KrigeST? Any ideas?


*SECOND EXPERIMENT Then I wondered what would happen if I just used distance as predictor so I coded instead: varPM10 <- variogramST(PM10~1,data=DFPM,tunit="days", tlags=0:0, assumeRegular=F,na.omit=T)

is this a reasonable way to try this? variogram and its model looks like:

image

using sepVgm <- vgmST("separable",space=vgm(1,"Per", 8, 700),time =vgm(200,"Exp", 15, 700), sill=100)

And the output really surprised me: image

why is this happening? I know the variogram modelling is poor, but even if that's true I understand the program should use the station data of the corresponding date so at least it should change in time.

My second question is: how would you fit the previous wave shaped variogram? I used ''show.vgms()'', and it appears that no model looks like that.


Thanks you in advance for any helpful advice, I apologize if I am asking something obvious.

Andrés Leslie