Объект не найден в R код

Question

Я пытаюсь выполнить следующий код, но в последних строках я нашел ошибку, потому что объект «HRdem» не найден (строка 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) 

Answer 1

Как лучше всего, как я могу сказать, что проблема начинается гораздо раньше:

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

В первой строке, я считаю, что x и y значение over заявлений перевернуты, поэтому она возвращает пустой столбец. Должно быть over(IDSTA.utm, grids). См. ?over в пакете sp для проверки того, что является.

Вторая строка не работает, потому что объект IDSTA.ov не имеет слота @data, потому что это не объект S4. Это просто обычная data.frame. over возвращает только data.frames или lists Возможно, вам понадобится заново создать его как объект SpatialGridDataFrame, с слотом data, который он создает по пути, прежде чем вы сможете использовать его так, как вы это делаете в остальной части программы. См. ?"SpatialGridDataFrame-class".

В результате этих проблем объект locs создан неправильно. Итак, когда вы вызываете функцию lm позже, она терпит неудачу при первом значении, которое бывает HRdem.

Надеюсь, это поможет.