This code creates the empirical semivariograms for three aerial images showing agricultural parcels. Want more info? Visit our www.gisandchips.org blog (http://www.gisandchips.org/2010/01/28/extraccion-de-caracteristicas-estructurales-de-una-imagen-ii-semi-variograma-de-una-imagen-usando-r-2/#more-1303).

empirical_semivariogram.R

# Instalo los packages necesarios: 
#-------------------------------------------------
 
install.packages("pixmap", dependencies= T)
 
install.packages("fields", dependencies= T)
 
# Cargo las librerias:
#-------------------------------------------------
library(pixmap)
 
library(fields)
 
# Hago una lista con las imagenes de las parcelas y las cargo todas:
#-------------------------------------------------
imag_dir <- list.files("C:\\ ... \\Articulos Gisandchips\\ parcelas_ppm\\", full.names=T)
 
parcela1 <- read.pnm(imag_dir[1])
 
parcela2 <- read.pnm(imag_dir[2])
 
parcela3 <- read.pnm(imag_dir[3])
 
# Aislamos una banda...
#-------------------------------------------------
matriz1<-parcela1@green
 
matriz2<-parcela2@green
 
matriz3<-parcela3@green
 
# Calculamos el semivariograma (si queremos podemos acabar antes cambiando el alcance de 40  a 10, por ejemplo, pero los resultados pueden cambiar y quedará menos bonito (-: )
#-------------------------------------------------
vgram1_40<-vgram.matrix( matriz1, R=40) # esto llevará un poco de tiempo
 
vgram2_40<-vgram.matrix( matriz2, R=40) # esto llevará un poco de tiempo
 
vgram3_40<-vgram.matrix( matriz3, R=40) # esto llevará un poco de tiempo
 
# Añadimos al Plot las matrices que acabamos de calcular
#-------------------------------------------------
plot.vgram.matrix(vgram1_40)# La matriz del variograma
 
plot.vgram.matrix(vgram2_40)# La matriz del variograma
 
plot.vgram.matrix(vgram3_40)# La matriz del variograma
 
# Creamos una curva que ajuste bien sobre la muestra y la añadimos al variograma
#-------------------------------------------------
polyfit1_40_20 <- lm(vgram1_40$vgram ~ poly(vgram1_40$d, 20));
 
plot(vgram1_40$d, vgram1_40$vgram)
 
lines(sort(vgram1_40$d), polyfit1_40_20$fit[order(vgram1_40$d)], col=2, lwd=4)
 
polyfit2_40_20 <- lm(vgram2_40$vgram ~ poly(vgram2_40$d, 20));
 
plot(vgram2_40$d, vgram2_40$vgram)
 
lines(sort(vgram2_40$d), polyfit2_40_20$fit[order(vgram2_40$d)], col=2, lwd=4)
 
polyfit3_40_20 <- lm(vgram3_40$vgram ~ poly(vgram3_40$d, 20));
 
plot(vgram3_40$d, vgram3_40$vgram)
 
lines(sort(vgram3_40$d), polyfit3_40_20$fit[order(vgram3_40$d)], col=2, lwd=4)