在 R 代码中找不到对象
Object Not found in R code
我正在尝试执行下一段代码，但在最后几行中我发现了错误，因为未找到对象 '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 (0C)", ylim=c(-12,28), main="GL001", col="grey")
scatter.smooth(KL003$cday, KL003$MDTEMP, xlab="Cumulative days", ylab="Mean daily temperature (0C)", ylim=c(-12,28), main="KL003", col="grey")
scatter.smooth(KL094$cday, KL094$MDTEMP, xlab="Cumulative days", ylab="Mean daily temperature (0C)", 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 (0C)", col="grey")
scatter.smooth(HRtemp2006locs$HRdsea, HRtemp2006locs$MDTEMP, xlab="Distance from the coast line (km)", ylab="Mean daily temperature (0C)", col="grey")
scatter.smooth(HRtemp2006locs$MODIS.LST, HRtemp2006locs$MDTEMP, xlab="MODIS LST (0C)", ylab="Mean daily temperature (0C)", 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)
据我所知，问题出现得更早：
IDSTA.ov <- over(grids, IDSTA.utm)
locs <- cbind(IDSTA.ov@data[c("HRdem", "HRdsea", "Lat", "Lon")], locs)
在第一行中，我相信 over 语句的 x 和 y 值是相反的，因此它 returns 是一个空列。应该是over(IDSTA.utm, grids)。请参阅 sp 包中的 ?over 以验证哪个是哪个。
第二行失败，因为 IDSTA.ov 对象没有 @data 槽，因为它不是 S4 对象。这只是一个普通的data.frame。 over 仅 returns data.frames 或 lists 您可能需要将其重新创建为 SpatialGridDataFrame 对象，并使用它沿在您可以按照您在程序其余部分中的方式使用它之前。有关详细信息，请参阅 ?"SpatialGridDataFrame-class"。
由于这些问题，locs 对象未正确创建。因此，当您稍后调用 lm 函数时，它会在第一个值处失败，该值恰好是 HRdem.
希望这对您有所帮助。
我正在尝试执行下一段代码，但在最后几行中我发现了错误，因为未找到对象 '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 (0C)", ylim=c(-12,28), main="GL001", col="grey")
scatter.smooth(KL003$cday, KL003$MDTEMP, xlab="Cumulative days", ylab="Mean daily temperature (0C)", ylim=c(-12,28), main="KL003", col="grey")
scatter.smooth(KL094$cday, KL094$MDTEMP, xlab="Cumulative days", ylab="Mean daily temperature (0C)", 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 (0C)", col="grey")
scatter.smooth(HRtemp2006locs$HRdsea, HRtemp2006locs$MDTEMP, xlab="Distance from the coast line (km)", ylab="Mean daily temperature (0C)", col="grey")
scatter.smooth(HRtemp2006locs$MODIS.LST, HRtemp2006locs$MDTEMP, xlab="MODIS LST (0C)", ylab="Mean daily temperature (0C)", 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)
据我所知，问题出现得更早：
IDSTA.ov <- over(grids, IDSTA.utm)
locs <- cbind(IDSTA.ov@data[c("HRdem", "HRdsea", "Lat", "Lon")], locs)
在第一行中，我相信 over 语句的 x 和 y 值是相反的，因此它 returns 是一个空列。应该是over(IDSTA.utm, grids)。请参阅 sp 包中的 ?over 以验证哪个是哪个。
第二行失败，因为 IDSTA.ov 对象没有 @data 槽，因为它不是 S4 对象。这只是一个普通的data.frame。 over 仅 returns data.frames 或 lists 您可能需要将其重新创建为 SpatialGridDataFrame 对象，并使用它沿在您可以按照您在程序其余部分中的方式使用它之前。有关详细信息，请参阅 ?"SpatialGridDataFrame-class"。
由于这些问题，locs 对象未正确创建。因此，当您稍后调用 lm 函数时，它会在第一个值处失败，该值恰好是 HRdem.
希望这对您有所帮助。