Search code examples
rsurfacetemperature

Object Not found in R code

I'm trying to execute the next code, but in the last lines I found an error because the object 'HRdem' is not found (line 161):

library(maptools)
library(gstat)
library(rgdal)
library(lattice)

# Download MODIS LST images:
download.file("http://spatial-analyst.net/book/sites/default/files/LST2006HR.zip", destfile=paste(getwd(), "LST2006HR.zip", sep="/"))
unzip(zipfile="LST2006HR.zip", exdir=getwd())
unlink("LST2006HR.zip")
download.file("http://spatial-analyst.net/book/sites/default/files/HRtemp2006.zip", destfile=paste(getwd(), "HRtemp2006.zip", sep="/"))
unzip(zipfile="HRtemp2006.zip", exdir=getwd())

HRtemp2006 <- read.delim("HRtemp2006.txt")
str(HRtemp2006) # Mean daily temperature for 365 days (2006) at 123 locations;

HRtemp2006$cday <- floor(unclass(as.POSIXct(HRtemp2006$DATE))/86400)
floor(unclass(as.POSIXct("2006-01-30"))/86400)[[1]]

IDSTA <- readShapePoints("IDSTA.shp", proj4string=CRS("+proj=longlat +datum=WGS84"))
IDSTA.utm <- spTransform(IDSTA, CRS("+proj=utm +zone=33 +ellps=WGS84
 + datum=WGS84 +units=m +no_defs"))
locs <- as.data.frame(IDSTA.utm)

names(locs) <- c("IDT_AK", "X", "Y")
str(locs)

dif.IDSTA <- merge(locs["IDT_AK"], data.frame(IDT_AK=levels(HRtemp2006$IDT_AK), sel=rep(1, length(levels(HRtemp2006$IDT_AK)))), by.x="IDT_AK", all.x=TRUE)

grids <- readGDAL("HRdem.asc")
names(grids@data)[1] <- "HRdem"
grids$HRdsea <- readGDAL("HRdsea.asc")$band1
proj4string(grids) <- IDSTA.utm@proj4string

grids.ll <- spTransform(grids[1], CRS("+proj=longlat +datum=WGS84"))
grids$Lat <- grids.ll@coords[,2]
grids$Lon <- grids.ll@coords[,1]
str(grids@data)

LST.listday <- dir(pattern=glob2rx("LST2006_**_**.LST_Day_1km.tif"))
LST.listnight <- dir(pattern=glob2rx("LST2006_**_**.LST_Night_1km.tif"))
for(i in 1:length(LST.listday)){
LSTname <- strsplit(LST.listday[i], ".LST_")[[1]][1]
tmp1 <- readGDAL(LST.listday[i])$band1
tmp2 <- readGDAL(LST.listnight[i])$band1
# convert to celsius:
tmp1 <- ifelse(tmp1<=7500, NA, tmp1*0.02-273.15)
tmp2 <- ifelse(tmp2<=7500, NA, tmp2*0.02-273.15)
grids@data[,LSTname] <- (tmp1+tmp2)/2

}
summary(grids$LST2006_05_17)

IDSTA.ov <- over(grids, IDSTA.utm)
locs <- cbind(IDSTA.ov@data[c("HRdem", "HRdsea", "Lat", "Lon")], locs)
str(locs)

HRtemp2006locs <- merge(HRtemp2006, locs, by.x="IDT_AK")
str(HRtemp2006locs)

LSTdate <- rep(NA, length(LST.listday))
for(i in 1:length(LST.listday)){
LSTdate[i] <- gsub("_", "-", strsplit(strsplit(LST.listday[i], ".LST_")[[1]][1], "LST")[[1]][2])
}
# cumulative days since 2006-01-01:
LSTcdate <- round((unclass(as.POSIXct(LSTdate))-unclass(as.POSIXct("2006-01-01")))/86400, 0) 
LSTcdate <- c(LSTcdate, 365)
LSTdate[1:5]; LSTcdate[1:5]

MODIStemp <- expand.grid(IDT_AK=levels(HRtemp2006$IDT_AK), DATE=levels(HRtemp2006$DATE), stringsAsFactors=TRUE)
MODIStemp$MODIS.LST <- rep(NA, length(MODIStemp[1]))
MODIStemp$MODIS.LST[1:(123*4)] <- rep(IDSTA.ov@data[!is.na(dif.IDSTA$sel), strsplit(LST.listday[i], ".LST_")[[1]][1]], 4)
for(i in 2:length(LST.listday)){
LSTname <- strsplit(LST.listday[i], ".LST_")[[1]][1]
d.days <- round((LSTcdate[i+1]-LSTcdate[i-1])/2, 0)
d.begin <- round((LSTcdate[i]-d.days/2)*123+1, 0)
d.end <- round((LSTcdate[i]+d.days/2)*123+1, 0)
MODIStemp$MODIS.LST[d.begin:d.end] <- rep(IDSTA.ov@data[!is.na(dif.IDSTA$sel),LSTname], d.days)
}
MODIStemp$MODIS.LST[(d.end+1):length(MODIStemp$MODIS.LST)] <- rep(IDSTA.ov@data[!is.na(dif.IDSTA$sel), strsplit(LST.listday[i], ".LST_")[[1]][1]], 2)

HRtemp2006locs$MODIS.LST <- MODIStemp$MODIS.LST[order(MODIStemp$IDT_AK)]
str(HRtemp2006locs)
tscale <- (((grids@bbox[1,"max"]-grids@bbox[1,"min"])+(grids@bbox[2,"max"]-grids@bbox[2,"min"]))/2)/(max(HRtemp2006locs$cday)-min(HRtemp2006locs$cday))
HRtemp2006locs$cdays <- tscale * HRtemp2006locs$cday
coordinates(HRtemp2006locs) <- c("X","Y","cdays")
proj4string(HRtemp2006locs) <- CRS(proj4string(grids))
HRtemp2006locs$c2days <- HRtemp2006locs@coords[,"cdays"]

MDTEMP.plt1 <- bubble(subset(HRtemp2006locs, HRtemp2006locs$cday==13150&!is.na(HRtemp2006locs$MDTEMP), select="MDTEMP"), fill=F, col="black", maxsize=2, key.entries=c(0,10,20,30), main="13150")
MDTEMP.plt2 <- bubble(subset(HRtemp2006locs, HRtemp2006locs$cday==13200&!is.na(HRtemp2006locs$MDTEMP), select="MDTEMP"), fill=F, col="black", maxsize=2, key.entries=c(0,10,20,30), main="13200")
MDTEMP.plt3 <- bubble(subset(HRtemp2006locs, HRtemp2006locs$cday==13250&!is.na(HRtemp2006locs$MDTEMP), select="MDTEMP"), fill=F, col="black", maxsize=2, key.entries=c(0,10,20,30), main="13250")
MDTEMP.plt4 <- bubble(subset(HRtemp2006locs, HRtemp2006locs$cday==13300&!is.na(HRtemp2006locs$MDTEMP), select="MDTEMP"), fill=F, col="black", maxsize=2, key.entries=c(0,10,20,30), main="13300")
print(MDTEMP.plt1, split=c(1,1,4,1), more=T)
print(MDTEMP.plt2, split=c(2,1,4,1), more=T)
print(MDTEMP.plt3, split=c(3,1,4,1), more=T)
print(MDTEMP.plt4, split=c(4,1,4,1), more=F)

GL001 <- subset(HRtemp2006locs@data, IDT_AK=="GL001", select=c("MDTEMP", "cday"))
KL003 <- subset(HRtemp2006locs@data, IDT_AK=="KL003", select=c("MDTEMP", "cday"))
KL094 <- subset(HRtemp2006locs@data, IDT_AK=="KL094", select=c("MDTEMP", "cday"))
par(mfrow=c(1,3))
scatter.smooth(GL001$cday, GL001$MDTEMP, xlab="Cumulative days", ylab="Mean daily temperature (\260C)", ylim=c(-12,28), main="GL001", col="grey")
scatter.smooth(KL003$cday, KL003$MDTEMP, xlab="Cumulative days", ylab="Mean daily temperature (\260C)", ylim=c(-12,28), main="KL003", col="grey")
scatter.smooth(KL094$cday, KL094$MDTEMP, xlab="Cumulative days", ylab="Mean daily temperature (\260C)", ylim=c(-12,28), main="KL094", col="grey")

par(mfrow=c(1,3))
scatter.smooth(HRtemp2006locs$HRdem, HRtemp2006locs$MDTEMP, xlab="DEM (m)", ylab="Mean daily temperature (\260C)", col="grey")
scatter.smooth(HRtemp2006locs$HRdsea, HRtemp2006locs$MDTEMP, xlab="Distance from the coast line (km)", ylab="Mean daily temperature (\260C)", col="grey")
scatter.smooth(HRtemp2006locs$MODIS.LST, HRtemp2006locs$MDTEMP, xlab="MODIS LST (\260C)", ylab="Mean daily temperature (\260C)", col="grey")

theta <- min(HRtemp2006locs$cday)
lm.HRtemp <- lm(MDTEMP~HRdem+HRdsea+Lat+Lon+cos((cday-theta)*pi/180)+MODIS.LST, data=HRtemp2006locs)
summary(lm.HRtemp)$adj.r.squared
hist(residuals(lm.HRtemp), col="grey", breaks=25)
plot(lm.HRtemp)
Solution

  • As best as I can tell, the problem starts much earlier:

    IDSTA.ov <- over(grids, IDSTA.utm)
locs <- cbind(IDSTA.ov@data[c("HRdem", "HRdsea", "Lat", "Lon")], locs)

    In the first line, I believe the x and y values of the over statement are inverted, so it returns a null column. Should be over(IDSTA.utm, grids). See ?over in the sp package to validate which is which.

    The second line fails because the IDSTA.ov object doesn't have an @data slot because it's not an S4 object. It's just a plain data.frame. over only returns data.frames or lists You may need to recreate it as a SpatialGridDataFrame object, with the data slot that it creates along the way before you can use it in the manner you have in the rest of the program. See ?"SpatialGridDataFrame-class" for details.

    As a result of these problems, the locs object is not created properly. So, when you call the lm function later on, it fails at the first value, which happens to be HRdem.

    Hope this helps.