Wheat Field Analysis

sampledata <-read.csv("Set4.196sample.csv")
wheatdata <- read.csv("Set4.196wheatyield.csv")

Data is sampled on a mostly regular grid with a few points missing in the bottom right. As such it alleviates challenges raised to random sampling by spatial data gathering processes. Interpolation would likely be valid.

coordinates(sampledata) <- c("Easting","Northing")
coordinates(wheatdata) <- c("Easting","Northing")

idwwheat <- idw(wheatdata$Yield~1, wheatdata, sampledata)

[inverse distance weighted interpolation]

spplot(idwwheat["var1.pred"])

The plot shows an inverse distance weight interpolation from the yield data to the sampling data. The plot suggests higher yields are observed at the bottom of the field.

head(idwwheat@data)

bdry.sf <- st_read("Set419697bdry.shp")

Reading layer `Set419697bdry' from data source `C:\Users\Localadmin\Documents\MATH944\Homework\Set419697bdry.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 592025 ymin: 4270352 xmax: 592470 ymax: 4271132
CRS:            NA

bdry.spdf <- as(bdry.sf, "Spatial")    

bbox(bdry.spdf)

      min     max
x  592025  592470
y 4270352 4271132

xseq <- seq(from= 592025, to = 592470, by= 5)
yseq <- seq(from = 4270352, to = 4271132, by= 5)
bdry.spdfdd <- expand.grid(xseq, yseq)
coordinates(bdry.spdfdd)<- c("Var1","Var2")

idwclay <- idw(sampledata$Clay~1, sampledata, bdry.spdfdd)

[inverse distance weighted interpolation]

spplot(idwclay["var1.pred"])

Gives an inverse distance weight interpolation from the sample data clay content onto a grid with distances of 5 meters between points. Plot suggests majority of clay content is at the top of the field.

clay.irr <- idwclay[bdry.spdf,]

spplot(clay.irr["var1.pred"])

Gives a subset of the previous plot including only the values included in the original irregular field.

data.var <- variogram(sampledata$Clay ~ 1, sampledata, cutoff= 600)
data.fit <- fit.variogram(data.var, model = vgm(1, "Gau", 700, 1))
plot(data.var,data.fit)

data.krig <- krige(sampledata$Clay~1, sampledata, bdry.spdfdd, model=data.fit)

[using ordinary kriging]

spplot(data.krig["var1.pred"])

The above is the krigging interpolation for clay based on a gausian variogram model. Results are consistent with the previous plots. Additionally a plot is included which displays the variogram model.

claymodel <-lm(Clay~ Northing, data=sampledata)
summary(claymodel)

Call:
lm(formula = Clay ~ Northing, data = sampledata)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.8470 -3.1530  0.5723  3.2655  6.0717 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.059e+04  8.248e+03  -10.98   <2e-16 ***
Northing     2.122e-02  1.931e-03   10.99   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.046 on 84 degrees of freedom
Multiple R-squared:  0.5897,    Adjusted R-squared:  0.5848 
F-statistic: 120.7 on 1 and 84 DF,  p-value: < 2.2e-16

plot(sampledata$Northing, sampledata$Clay)
abline(claymodel)

sampledataresid <-add_residuals(sampledata, claymodel, var ="resid")

I detrend the data with a SLR relating Clay to the northing coordinates.

resid.var <-variogram(sampledataresid$resid~1, sampledataresid, cutoff= 600)
resid.fit <- fit.variogram(resid.var, model = vgm(1, "Gau", 700, 1))

resid.krig <- krige(sampledataresid$resid~1, sampledataresid, bdry.spdfdd, model=resid.fit)

[using ordinary kriging]

plot(resid.var,resid.fit)
resid.krig$var2.pred <- coef(claymodel)[1] + coef(claymodel)[2]*bdry.spdfdd@coords[,2] + resid.krig$var1.pred
par(mfrow = c(1,2))

spplot(resid.krig["var2.pred"])

spplot(resid.krig["var1.pred"])

head(bdry.spdfdd)

class       : SpatialPoints 
features    : 1 
extent      : 592025, 592025, 4270352, 4270352  (xmin, xmax, ymin, ymax)
crs         : NA 

Krigging interpolation is done on the residuals and the model is retrended. Plot is provided for both residuals and retrended.

par(mfrow = c(1,2))

spplot(resid.krig["var2.pred"])

spplot(data.krig["var1.pred"])

Side by side comparison for trended and detrended data. Results are similar, But detrended more accurately reflects the concentration of clay content.