je souhaite faire un Krigeage ordinaire sous R.
j'ai créer ma grille à partir d'un vecteur du contour de ma région.
Code : Tout sélectionner
#import area border from ESRI shape file
border<-readOGR(dsn="D:/00-New architecture/2-RELEPA/9-Analyses/Krigeage","border")
proj4string(border)=CRS("+init=epsg:2154")
str(border)
#Let's first create a prediction grid for the interpolation, starting from the shape file
vals <- border@bbox
deltaLong <- as.integer((vals[1,2] - vals[1,1]) + 1.5)
deltaLat <- as.integer((vals[2,2] - vals[2,1]) + 1.5)
gridRes <-0.01 #change this value to change the grid size (in metres)
gridSizeX <- deltaLong / gridRes
gridSizeY <- deltaLat / gridRes
grd <- GridTopology(vals[,1],c(gridRes,gridRes),c(gridSizeX,gridSizeY))
pts <- SpatialPoints(coordinates(grd))
pts1 <- SpatialPixelsDataFrame(as.data.frame(pts), data=as.data.frame(rep(1,nrow(as.data.frame(pts)))))
Overlay=overlay(pts1,border)
pts1$border=Overlay
nona<-na.exclude(as.data.frame(pts1))
coordinates(nona)<-~x+y
gridded(nona) <- TRUE
proj4string(nona)=CRS("+init=epsg:2154") #remeber to set the coordinate system also for the prediction grid
writeAsciiGrid(nona,"prediction_grid.asc")
org_mat<-read.asciigrid("prediction_grid.asc")
proj4string(org_mat)=CRS("+init=epsg:2154")
jusque là ça va, qand j'affiche le plot ça correspond.
Ensuite j'ai fait mon variogramme de la manière suivante :
Code : Tout sélectionner
#Fitting the variogram
#first plot the variogram
vgm<-variogram(IND~1,data)
vgm
plot(vgm,type="o",pch=3)
#then, create a variogram model which is better(2)
mod2<-vgm(psill=var(data$IND),model="Exp",range=90000,nugget=0) #From Spatial-analyst.net
fit2<-fit.variogram(vgm,mod2)
plot(vgm,fit2,main="Model Exp",type="p",pch=10)## ok
pas d'erreur dans R jusque là.
par contre la suite me pose problème lorsque que je veux kriger mes résultats ma carte est unicolor comme si les prédiction étaient toutes identiques...
quand je lance le krigeage avec cette manip :
Code : Tout sélectionner
map<-krige(IND~1,data,nona,model=fit2)
spplot(map,"var1.pred",col.regions=heat.colors(50))
aucune erreur n'apparait sur la console R, mais c'est là que j'obtiens une carte unicolor.
est ce que quelqu'un sait d'où cela peut provenir ? (si besoin d'info supplémentaire sur la forme de mes données n’hésitez pas.
mes valeurs predites sont toutes identiques, et je ne vois pas dans mon script ce qui donne ce résultats.
Merci d'avance pour vos réponses et votre interet porter à ma question.
je n'ai pas réussi à trouver quelqu'un ayant eu le même problème.
Et meilleurs voeux à tous !