Geostatistik - regression Kriging

[DAFZ]

Liebes Forum,

ich brauch bitte eure Hilfe. Ich sitze nun schon seit drei Wochen an einem eigentlichen einfachen Problem und bekomme es nicht geloest. Ich wuerde gerne die Temperatur gegen die Erhebung auf ein Wassereinzugsgebiet interpolieren, bekomme aber immer wieder die gleichen Fehler. Ich weiss mir nicht mehr zu helfen. Ich bin ein ziemlicher Anfaenger in R und habe schon alles mir bekannte versucht und natuerlich auch das Internet zu Rate gezogen, aber ich kann das Problem nicht loesen.
Ich arbeite hauptsaechlich mit dem Paket "gstat". Ich habe 15 Messpunkte innerhalb und ausserhalb des Einzuggebietes und moechte ausgehend von diesen Messpunkten interpolieren. Das Einzugsgebiet habe ich ebenfalls als Dataframe in Form von Punktkoordinaten (315226 Reihen) mit der Erhebung als zusaetzlicher Information.

Mein Code ist wie folgt:

Code: Alles auswählen

library(sp)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)
library(gstat)

d_t<-data_temperature
coordinates(d_t)= ~LONG+LAT

jh_sp<- jh_coord
coordinates(jh_sp)=~X+Y

ext_jh<-extent(jh_sp)

jh_ras<- raster(ext_jh)

jh_raster<-rasterize(jh_sp, jh_ras,jh_sp$GRID_CODE, fun='last', background=NA,
          mask=FALSE, update=FALSE, updateValue='all', na.rm=TRUE)

proj4string(jh_raster)<-CRS("+init=epsg:4326")

jh_grid <- rasterToPoints(jh_raster, spatial = TRUE)

gridded(jh_grid)<- TRUE

proj4string(d_t) <- CRS("+init=epsg:4326")

stat_prj_t <- spTransform(d_t, CRSobj = CRS("+init=epsg:4326"))

lzn.vgm = variogram(log(X1)~1, stat_prj_t)

lzn.fit = fit.variogram(lzn.vgm, model = vgm(1, "Sph", 4000, 1))

plot(lzn.vgm, lzn.fit)

lzn.kriged = krige(stat_prj_t$X1~stat_prj_t$ELEVATION , newdata=stat_prj_t, locations=jh_grid, model = lzn.fit)

Anstelle habe ich auch folgendes ausprobiert:

Code: Alles auswählen

gjh <- gstat(id = "stat_prj_t", formula = stat_prj_t$X1~stat_prj_t$ELEVATION,
                data = stat_prj_t, maxdist=100, nmin=10, force=TRUE, model=lzn.fit)

predgjh<- predict(gjh, jh_grid)

Wenn ich:

Code: Alles auswählen

lzn.kriged = krige(stat_prj_t$X1~stat_prj_t$ELEVATION , newdata=stat_prj_t, locations=jh_grid, model = lzn.fit)

ausfuehre, bekomme ich folgenden Fehler:

Code: Alles auswählen

Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
      dimensions do not match: locations 154 and data 15

Wenn ich:

Code: Alles auswählen

 gjh <- gstat(id = "stat_prj_t", formula = stat_prj_t$X1~stat_prj_t$ELEVATION,
                    data = stat_prj_t, maxdist=100, nmin=10, force=TRUE, model=lzn.fit)
    predgjh<- predict(gjh, jh_grid)

ausfuehre, bekomme ich diesen Fehler:

Code: Alles auswählen

 Error in gstat.formula.predict(d$formula, newdata, na.action = na.action,  : 
      NROW(locs) != NROW(X): this should not occur

Was bedeutet das? Wo liegt der Fehler in meinen Daten? Wie schaffe ich es, dass sie die gleichen Dimensionen besitzen und das Kriging problemlos verlaeuft?

[EDi]

Bitte ein reproduzierbares Beispiel posten, sonst muss man raten was passiert...

Ich würde raten:

Code: Alles auswählen

lzn.kriged = krige(X1 ~ ELEVATION , locations = stat_prj_t, newdata = jh_grid, model = lzn.fit)

Hier ein reproduzierbares Beispiel:

Code: Alles auswählen

library(sp)
library(gstat)
data("meuse")
# convert to spatialPointsDataframe
coordinates(meuse) <- ~ x + y
plot(meuse)

# fit a variogram
lzn.vgm <- variogram(log(zinc)~1, meuse) 
lzn.fit <- fit.variogram(lzn.vgm, model=vgm(1, "Sph", 900, 1)) # fit model
plot(lzn.vgm, lzn.fit)

# grid where to interpolate
data("meuse.grid")
coordinates(meuse.grid) <- ~ x + y # step 3 above
plot(meuse.grid)

# krige
lzn.kriged <- krige(log(zinc) ~ 1, locations = meuse, newdata = meuse.grid, model = lzn.fit)
spplot(lzn.kriged["var1.pred"])

[EDi]

Generell kann man sich merken, dass '$' in formulas in den allermeisten Fällen keine gute Idee sind (wie immer gibt es Ausnahmen).

[DAFZ]

Danke EDI

Das mit dem Beispiel hatte ich vergessen. Tut mir Leid. Das Beispiel ist aus dem gstat Paket oder?

Wenn ich dein Ratschlag befolge, gibt er mir an ELEVATION nicht zu finden. ELEVATION wird aber genauso unter View(stat_prj_t) mit aufgefuehrt.
Kannst du mir erklaeren worauf sich "new data =" und "locations=" beziehen?
Ich habe immer gedacht, dass locations angibt, auf welche Flaeche interpoliert werden soll.

[EDi]

Steht doch in der Hilfe... new_data ist die Fläche auf interpoliert werden soll...

Wenn ich dein Ratschlag befolge, gibt er mir an ELEVATION nicht zu finden. ELEVATION wird aber genauso unter View(stat_prj_t) mit aufgefuehrt.

Ein reproduzierbares Beispiel fehlt immernoch, deshalb kann ich hier nichts mehr zu sagen.

Das Beispiel ist aus dem gstat Paket oder?

Folgt grob

Code: Alles auswählen

example(krige)

[DAFZ]

Vielen Dank! Das hat mir schon sehr geholfen.

Ich habe das Problem letztendlich geloest, in dem ich kein Grid gebaut habe, sondern einfach auf die Punktdaten meines Einzuggebietes interpoliert habe.