[R-sig-Geo] climate data interpolation

[Mintewab Bezabih]

Dear the R-sig-geo community,


I have been trying to interpolate rainfall and temperature data for 8
stations and that needs to be interpolated at a farm level so that I
have farm level temperature and rainfall observations.

I have attempted handling the problem with the code below and while I
succedded in getting results, my values look very wrong.

I am now revisiting what I did and I realized that perhaps it is 3d
kriging in gstat that I have to do. Since my grasp of R is not that
great, I just thought of asking you for tips on how I do a correct
interpolation.

Below is a reproducable example of what I did

many thanks
minti


#library(mgcv)
library(fields)

x <-1:20
y<- runif(20)
z<- c(11, 15, 17, 2, 18, 6, 7, NA, 12, 10,21, 25, 27, 12, 28, 16,
17, NA, 12, 10)
d<- c(11, 15, 17, 2, 18, 6, 7, NA, 12, 10,21, 25, 27, 12, 28, 16,
17, NA, 12, 10)

z <- dat$sum360
x <- dat$latitud1
y <- dat$longtud1
d <- dat$seqnum
d1 <- dat$seqnum
d2 <- dat$seqnum
d[ is.na(x) ] <- NA
d[ is.na(y) ] <- NA
d[ is.na(z) ] <- NA


F<-data.frame(x, y, z,d, d1,d2)
mydataset<-data.frame(x, y, z,d)
airquality<-mydataset[!is.na(x) & !is.na(y)& !is.na(z)& !is.na(d), ]
airquality[complete.cases(mydataset),]


x<- airquality[, c(1)]
y<- airquality[, c(2)]
z<- airquality[, c(3)]
d1<- airquality[, c(4)]

#xmin <- min(x)
#xmax <- max(x)
#ymin <- min(y)
#ymax <- max(y)


xmin <- min(x) - 1/2
xmax <- max(x) + 1/2
ymin <- min(y) - 1/2
ymax <- max(y) + 1/2


F1<-data.frame(x, y,z)

tpsfit <- Tps(cbind(x, y), z, scale.type="unscaled")
#tpsfit2 <- Tps(cbind(x, y), z, lambda=20*tpsfit$lambda,scale.type="unscaled")


surface (tpsfit)
out.p<- predict(tpsfit)

library(raster)
#setwd("f:/")

#setwd("f:/climate change data/climateChange paper/thin plate/raster2/raster")
#r2 <- raster("current_bio_1_1.asc")
#res(r2) <- 0.00008
#r2 <- raster("amhara_coordinates GCS_WGS 1984.asc")
r <- raster(as.matrix(F1))
extent(r) <- extent(xmin, xmax, ymin, ymax)
#r2 <- raster("raster_asc.txt")
#r2 <- expand(r2, extent(r2)+0.175)
#newExtent <- extent(c(60, 100, 0, 40))
#newExtent <- extent(c(30, 45, 10, 40))
#newExtent <- extent(c(0, 84, 0, 84))
r2 <- crop(r, extent(r))

plot(r2)

#ra <- aggregate(r, 10)
#xy <- data.frame(xyFromCell(ra, 1:ncell(ra)))
#v <- getValues(ra)
#tps <- Tps(xy, v)
#m <- raster(r2)
p <- interpolate(r2, tpsfit)
#p <- mask(p, r2)
plot(p)
#se <- interpolate(p, tpsfit, fun=predict.se)
#se <- mask(se, r3)
#plot(se)

#xy <- data.frame(x,y,z)
require(sp)
head(as.data.frame(as(p, "SpatialGridDataFrame")))
result1 <- as.data.frame(as(p, "SpatialGridDataFrame"))
result1 <- data.frame(result1, d1)
#names(result1)[names(result1)=="d1"]<-"d1"
#result2<-data.frame(d1=row.names(result1),result1[, c("layer","s1","s2")])
result2<-merge(F, result1,
by=c("d1"),all.x=TRUE,all.y=TRUE,all.z=TRUE,all.layer=TRUE)
#result2<-merge(F, result1, by=c("d1"))

#result3<-data.frame(result2[ , c("x","y","z", "layer")])

#result3<-result2[!duplicated(result2[,c("x","y","z","d1")]),]

result3<-result2[!duplicated(result2[,c("d1")]),]

#result3<-data.frame(result3[ , c("x","y","z", "layer")])

#result2<-merge(mydataset, result1, by=c("x","y"),all.x=TRUE,all.y=TRUE)
write.table(result3, file = "f:/climate change data/climateChange
paper/thin plate/interpolated_precip.csv", sep = ",", col.names = NA)