La interpolación de Thin Plate Spline corresponde a las técnicas de interpolación n-dimensional a partir de datos dispersos. La función de interpolación correspondiente es una función continuamente diferenciable (función C1) que tiene la energía (norma L2) de su segundo derivado mínimo.
El residuo de interpolación con respecto a esta interpolación puede ser 0 ó, lo que es más interesante, se puede hacer una compensación de regularidad entre el error de interpolación y la energía de la curvatura de la función.
Wahba (1990) propone un método basado en la validación cruzada para determinar la compensación óptima.
- Thin Plate Spline adapta los datos espaciados irregularmente a una superficie.
- El parámetro de suavizado se elige mediante validación cruzada generalizada.
- El modelo asumido es el aditivo Y = f (X) + e donde f (X) es una superficie dimensional.
- Este es un caso especial de la estimación espacial.
Thin Plate Spline es el resultado de minimizar la suma residual de cuadrados sujetos a una restricción de que la función tiene un cierto nivel de suavidad (o penalización de rugosidad). En donde la rugosidad se cuantifica por la integral de los derivados de orden m cuadrado al cuadrado.
Para una dimensión y m = 2, la penalización de rugosidad es el cuadrado integrado de la segunda derivada de la función. Para dos dimensiones, la penalización por rugosidad es la integral de
(Dxx (f))^2 + 2 (Dxy (f)) ^2 + (Dyy (f)) ^2
Además de controlar el orden de las derivadas, el valor de m también determina el polinomio base que se ajusta a los datos. El grado de este polinomio es (m-1).
El parámetro de suavizado controla la cantidad de suavización de los datos. En la forma habitual, esto se denota con lambda, el multiplicador de Lagrange del problema de minimización. Aunque esta es una escala incómoda, lambda = 0 no corresponde con restricciones de suavidad y los datos se interpolan. lambda = infinito corresponde a simplemente ajustar el modelo de base polinomial por mínimos cuadrados ordinarios.
Este estimador se implementa pasando la función de covarianza generalizada correcta basada en funciones de base radial a la función más general Krig. Una ventaja de esta implementación es que una vez que se crea un objeto Tps / Krig, el estimador se puede encontrar rápidamente para otros datos y parámetros de suavizado siempre que las ubicaciones permanezcan sin cambios. Esto hace que la simulación dentro de R sea eficiente. Tps no admite actualmente el argumento de los nudos donde se puede usar un conjunto reducido de funciones básicas. Esto es principalmente para simplificar y una buena alternativa al usar nudos sería usar una covarianza válida de la familia Matern y un parámetro de gran rango.
Código utilizado en R
##################################################################
#
# Walter Arnoldo Bardales Espinoza (bardaleswa@gmail.com)
# julio de 2018
#
##################################################################
# Cargar librerias
library(fields)
library(raster)
library(rgdal)
library(lattice)
library(stringr)
# Cargar modelo de elevación
dem = raster("D:/Documents/GIS DataBase/DEM_CA/geotif/dem_90_gtm.tif")
# Extraer las coordenadas del DEM
coords = xyFromCell(dem,1:ncell(dem))
# Extraer los factores de las coordenadas x e y
coord_x = as.numeric(levels(factor(coords[,"x"])))
coord_y = as.numeric(levels(factor(coords[,"y"])))
# Convertir el dem en matriz y transponerlo
z = t(as.matrix(dem))
# Corregir la matriz de elevacion debido a que esta rotada 180 grados
coord_z = matrix(, nrow=length(coord_x), ncol=length(coord_y))
for(i in 1:length(coord_y)){
j = length(coord_y)+1-i
coord_z[,j] = z[,i]
}
#crear lista de elevación a partir del dem
elevacion = list(x=coord_x,y=coord_y,z=coord_z)
# Evaluar la lista
str(elevacion)
# Cargar los datos de lluvia
datospp = read.csv("D:/Documents/Lluvia/Lluviaxy_1980_2016.csv",sep=",")
datos = datospp
# Convertir en puntos geograficos y asignarle las coordenadas GTM
coordinates(datos) = ~x + y
gtm = CRS("+proj=tmerc +lat_0=0 +lon_0=-90.5 +y_0=0 +x_0=500000 +k=0.9998 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj4string(datos) = gtm
# Extraer los valores de elevacion de las estaciones a partir del DEM
altura = extract(dem,datos)
# Crear lista con los datos x, y, z y variable
xydatos = data.frame(lon=datospp$x,lat=datospp$y)
zdatos = altura
# Antes de continuar verificar que no hallan datos faltantes en las series
# Nombrar columnas
mes = c(0,0,0,1:13)
meses = c("","","","enero", "febrero", "marzo", "abril", "mayo", "junio", "julio", "agosto", "septiembre", "octubre", "noviembre", "diciembre", "Anual")
for(k in 4:16){
ydatos = datospp[,k]
precip = list(x=xydatos,elev=zdatos,y=ydatos)
# Graficar las estaciones y valores de lluvia
quilt.plot(precip$x,precip$y)
# Crear el modelo de lluvia elevacion
fit = Tps(precip$x, precip$y, Z=precip$elev)
# Resumen del modelo lluvia elevacion con TPS
summary(fit)
# Graficar los resultados
set.panel(2,2)
plot(fit)
# Intrpolar superficie del modelo creado
set.panel()
# Crear la lista de coordenadas a interpolar
grid.list = list(x=elevacion$x, y=elevacion$y)
#Interpolacion del modelo usando de covariable la elevación
fit.full = predictSurface(fit, grid.list, extrap=TRUE, ZGrid = elevacion)
# Transponer la matriz de resultados
pp_z1 = t(fit.full$z)
# Crear el raster interpolado
ras = raster(pp_z1, xmn=extent(dem)[1],xmx=extent(dem)[2],ymn=extent(dem)[3],ymx=extent(dem)[4])
ras = flip(ras,direction="y")
crs(ras) = gtm
ras[ras<0 0="" p="">
# Guardaer el raster
destino5 = paste("D:/Documents/Consultorias_2018/balance_hidrico/Raster/lluvia/prueb/prec_",str_pad(mes[k],width =2, pad="0"),".tif", sep ="")
writeRaster(ras,filename=destino5, overwrite=TRUE)
}
Resultados de la interpolación
