HRtemp08: Daily temperatures for Croatia for year 2008

Description

The daily measurements of temperature (thermometers) for year 2008 kindly contributed by the Croatian National Meteorological Service. HRtemp08 contains 56,608 measurements of temperature (159 stations by 365 days).

Usage

data(HRtemp08)

Arguments

Format

The HRtemp08 data frames contain the following columns:

NAME: name of the meteorological station
Lon: a numeric vector; x-coordiante / longitude in the WGS84 system
Lat: a numeric vector; y-coordinate / latitude in the WGS84 system
DATE: 'Date' class vector
TEMP: daily temperature measurements in degree C

References

Hengl, T., Heuvelink, G.B.M., Percec Tadic, M., Pebesma, E., (2011) Spatio-temporal prediction of daily temperatures using time-series of MODIS LST images. Theoretical and Applied Climatology, 107(1-2): 265-277. 10.1007/s00704-011-0464-2
AGGM book datasets (http://spatial-analyst.net/book/HRclim2008)

Examples

Run this code

# NOT RUN {
data(HRtemp08)

# }
# NOT RUN {
## examples from: http://dx.doi.org/10.1007/s00704-011-0464-2 
library(spacetime)
library(gstat)
library(sp)
sp <- SpatialPoints(HRtemp08[,c("Lon","Lat")])
proj4string(sp) <- CRS("+proj=longlat +datum=WGS84")
HRtemp08.st <- STIDF(sp, time = HRtemp08$DATE-.5, 
     data = HRtemp08[,c("NAME","TEMP")], 
     endTime = as.POSIXct(HRtemp08$DATE+.5))
## Country borders:
con0 <- url("http://www.gadm.org/data/rda/HRV_adm1.RData")
load(con0)
stplot(HRtemp08.st[,"2008-07-02::2008-07-03","TEMP"], 
   na.rm=TRUE, col.regions=SAGA_pal[[1]],
   sp.layout=list("sp.polygons", gadm))

## Load covariates:
con <- url("http://plotkml.r-forge.r-project.org/HRgrid1km.rda")
load(con)
str(HRgrid1km)
sel.s <- c("HRdem","HRdsea","HRtwi","Lat","Lon")
## Prepare static covariates:
begin <- as.Date("2008-01-01")  
endTime <- as.POSIXct(as.Date("2008-12-31"))
sp.grid <- as(HRgrid1km, "SpatialPixels")
HRgrid1km.st0 <- STFDF(sp.grid, time=begin, 
    data=HRgrid1km@data[,sel.s], endTime=endTime)
## Prepare dynamic covariates:
sel.d <- which(!names(HRgrid1km) %in% sel.s)
dates <- sapply(names(HRgrid1km)[sel.d], 
     function(x){strsplit(x, "LST")[[1]][2]}
)
dates <- as.Date(dates, format="%Y_%m_%d")
## Sort values of MODIS LST bands:
m <- data.frame(MODIS.LST = as.vector(unlist(HRgrid1km@data[,sel.d])))
## >10M values!
## Create an object of type STFDF:
HRgrid1km.stD <- STFDF(sp.grid, time=dates-4, data=m, 
     endTime=as.POSIXct(dates+4))

## Overlay in space and time:
HRtemp08.stxy <- spTransform(HRtemp08.st, CRS(proj4string(HRgrid1km)))
ov.s <- over(HRtemp08.stxy, HRgrid1km.st0)
ov.d <- over(HRtemp08.stxy, HRgrid1km.stD)
## Prepare the regression matrix:
regm <- do.call(cbind, list(HRtemp08.stxy@data, ov.s, ov.d))
## Estimate cumulative days:
regm$cday <- floor(unclass(HRtemp08.stxy@endTime)/86400-.5)
str(regm)
## Plot a single station:
scatter.smooth(regm$cday[regm$NAME=="Zavi<c5><be>an"], 
    regm$TEMP[regm$NAME=="Zavi<c5><be>an"], 
    xlab="Cumulative days", 
    ylab="Mean daily temperature (\260C)", 
    ylim=c(-12,28), main="GL039 (Zavi\236an)", 
    col="grey")
## Run PCA so we can filter missing pixels in the MODIS images:
pca <- prcomp(~HRdem+HRdsea+Lat+Lon+HRtwi+MODIS.LST, 
    data=regm, scale.=TRUE)
selc <- c("TEMP","Lon","Lat","cday")
regm.pca <- cbind( regm[-pca$na.action, selc], 
    as.data.frame(pca$x))
## Fit a spatio-temporal regression model:
theta <- min(regm.pca$cday)
lm.HRtemp08 <- lm(TEMP~PC1+PC2+PC3+PC4+PC5+PC6 
     +cos((cday-theta)*pi/180), data=regm.pca)
summary(lm.HRtemp08)

## Prediction locations -> focus on Istria:
data(LST)
gridded(LST) <- ~lon+lat
proj4string(LST) <- CRS("+proj=longlat +datum=WGS84")
LST.xy <- reproject(LST[1], proj4string(HRgrid1km))
LST.xy <- as(LST.xy, "SpatialPixels")
## targeted dates:
t.dates <- as.Date(c("2008-02-01","2008-05-01","2008-08-01"), 
     format="%Y-%m-%d")
LST.st <- STF(geometry(LST.xy), time=t.dates)
## get values of covariates:
ov.s.IS <- over(LST.st, HRgrid1km.st0)
ov.d.IS <- over(LST.st, HRgrid1km.stD)
LST.stdf <- STFDF(geometry(LST.xy), time=t.dates, 
    data=cbind(ov.s.IS, ov.d.IS))
## predict Principal Components:
LST.pca <- as.data.frame(predict(pca, LST.stdf@data))
LST.stdf@data[,paste0("PC",1:6)] <- LST.pca
cday.l <- as.vector(sapply(
    floor(unclass(LST.stdf@endTime)/86400-.5), 
    rep, nrow(LST.xy@coords)))
LST.stdf@data[,"cday"] <- cday.l
stplot(LST.stdf[,,"PC1"], col.regions=SAGA_pal[[1]])
stplot(LST.stdf[,,"PC2"], col.regions=SAGA_pal[[1]])

## Predict spatio-temporal regression:
LST.stdf@data[,"TEMP.reg"] <- predict(lm.HRtemp08, 
     newdata=LST.stdf@data)
## Plot predictions:
gadm.ll <- as(spTransform(gadm, 
     CRS(proj4string(HRgrid1km))), "SpatialLines")
stplot(LST.stdf[,,"TEMP.reg"], col.regions=SAGA_pal[[1]],
  sp.layout=list( list("sp.lines", gadm.ll),
   list("sp.points", HRtemp08.stxy, col="black", pch=19) )
  )
# }