Skip to content

[R-es] Duda interpolación (package ' gstat ')

6 messages · Marcos Bermejo, Freddy López, rubenfcasal

#
Hola,

  # Hacemos el KED. Ver función "krige()":
        KED.rad <- krige(
          formula=pluvPcp~layer,                  # covariable -> radar
          locations=lluvia.rad.pluv.spdf,
          newdata=radarGrid,                      # podría ser cualquier objeto Spatial
          model=v.fit,                            # modelo de semivariograma. 
          maxdist=Inf
        )

Esta es la función que me interpola los datos de lluvia. El error que me da es:

"solve.c", line 88: singular matrix in function Usolve()

"lufactor.c", line 208: singular matrix in function m_inverse()
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
  m_inverse
In addition: Warning message:
In fit.variogram(vg.aux, model = vgm(psill = 0.1, model = "Gau",  :
  Warning: singular model in variogram fit


Mi función del variograma es :
v.fit <- fit.variogram(vg.aux, model=vgm(psill=0.15, model='Gau', range=5000,
                                               nugget=0.05))

¿Alguien me podría ayudar?

Gracias de antemano.

Un saludo,
2 days later
#
Hola Marcos,

     Parece que el problema es con el ajuste del variograma (sale 
plano?), sin más información no se exactamente que puede estar pasando...

     Si me envías el código completo y los datos lo miro con más detalle 
(e incluso te doy una alternativa no paramétrica con el paquete npsp).

     Un saludo, Rubén.



El 04/08/2015 a las 11:24, Marcos Bermejo escribió:

  
  
#
Sale plano sí. 

Ya se que sin tener los datos y el código es un poco difícil, pero es que mis datos ocupan mucho, es imposible.
Seguiré mirando por internet. 

Muchas gracias Rubén.

Un saludo,

  
  
#
Hola Marcos,

¿El problema persiste si pruebas con un subconjunto de los datos?

Saludos.

2015-08-06 12:34 GMT-03:00 Marcos Bermejo <markbermejo90 en hotmail.com>:

  
    
#
Hola de nuevo,

     Un par de detalles técnicos...

     Antes de nada comentar que el paquete gstat no es 
computacionalmente muy eficiente calculando las predicciones kriging 
(calcula el estimador mcg de la tendencia utilizando la expresión 
explícita, etc...), entre otras cosas requiere la factorización de la 
matriz de varianzas covarianzas y pueden aparecer problemas numéricos.

     Es bien conocido que con el modelo gaussiano de variograma pueden 
aparecer estas inestabilidades (puede ser muy plano en saltos pequeños y 
como consecuencia  la matriz de covarianzas es semidefinida positiva 
pero no 'estrictamente' definida positiva - en esto influye el redondeo...).

     Mi recomendación sería que probases con otro modelo de variograma y 
que compartas el gráfico del ajuste...

     Un saludo, Rubén.

P.D. Cuidado también con el sesgo en la estimación del variograma a 
partir de los residuos (e.g. Fernandez-Casal R. and Francisco-Fernandez 
M. (2014) Nonparametric bias-corrected variogram estimation under 
non-constant trend, Stoch. Environ. Res. Ris. Assess, 28, 1247-1259), 
aunque si tu objetivo final es la predicción no te preocupes demasiado 
(no deberías fiarte de las varianzas kriging)...


El 06/08/2015 a las 17:40, Freddy Omar López Quintero escribió:
1 day later
#
Muchas gracias por todo.

Esta es la parte del código del variograma: 
vg.aux <- variogram(radPcp~1, radar.spdf, cutoff=10000)
vg <- as.data.frame(matrix(c(vg.aux$dist, vg.aux$gamma), nrow=15, ncol=2))
# El primer argumento de "fit.variogram()" es lo que hemos obtenido de "variogram()".
v.fit <- fit.variogram(vg.aux, model=vgm(psill=0.10, model='Gau', range=5000,
                                               nugget=0.05))

He probado con otros modelos: 'Exp', 'Mat' y 'Sph', y me sigue dando el mismo error que os mandé.

Os pongo la función "krige()" otra vez:
KED.rad <- krige(
          formula=pluvPcp~layer,                  # covariable -> radar
          locations=lluvia.rad.pluv.spdf,
          newdata=radarGrid,                      # podría ser cualquier objeto Spatial
          model=v.fit,                            # modelo de semivariograma. 
          maxdist=Inf
        )
        

La cosa es que tengo datos de lluvia de 5 pluviómetros (variable 'pluvPcp') situados en 5 puntos diferentes de una cuenca, y quiero interpolarpolarlos para que me estime la lluvia en toda la cuenca usando como covariable los datos de un radar meteorológico en esos mismos 5 puntos de los pluviómetros (variable 'layer'). 
El variograma lo obtengo a partir de los datos del radar en toda la cuenca (variable 'radPcp' -> unos 13.000 puntos), ya que no puedo obtenerlo a partir de los datos de lluvia de los pluviómetros, porque con 5 puntos no puedo obtener un variograma fiable.
¿El error que me daba R, puede ser por esto? ¿Cómo lo veis?

Lo de la solución no paramétrica no lo entiendo muy bien y no sé si podré implementarlo en mi caso.

Un saludo,