Mapping spatial data

#install.packages("ggmap")
suppressPackageStartupMessages({
    library(ggmap)
    library(dplyr)
    library(scales)
    library(magrittr)
})

Load map of WWPS

WWPS_map <- get_map(location=c(lon=-94.2348,lat=36.0661),zoom=19,maptype = "satellite",source="google")

Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=36.0661,-94.2348&zoom=19&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false

ggmap(WWPS_map)

Read in data table with sample positions

sample_locations <- read.csv("C:/Users/faysmith/Desktop/Export_Output.csv")

Plot locations on map

sam_map_gwc <- ggmap(WWPS_map) + geom_point(data = sample_locations, aes(lon,lat,color=gwc), size=2,alpha=0.7)+labs(x="longitude", y="latitude", title = "WWPS Soil Gravimetric Water Content 0-10cm", color="GWC (g H20/g soil)")
sam_map_gwc

sam_map_ph <- ggmap(WWPS_map) + geom_point(data = sample_locations, aes(lon,lat,color=ph1), size=2,alpha=0.7)+labs(x="longitude", y="latitude", title = "WWPS Soil pH 0-10cm", color="pH")
sam_map_ph

sam_map_ec <- ggmap(WWPS_map) + geom_point(data = sample_locations, aes(lon,lat,color=ec1), size=2,alpha=0.7)+labs(x="longitude", y="latitude", title = "WWPS Soil EC 0-10cm", color="EC (uS/cm)")
sam_map_ec

sam_map_clay <- ggmap(WWPS_map) + geom_point(data = sample_locations, aes(lon,lat,color=per_clay), size=2,alpha=0.7)+labs(x="longitude", y="latitude", title = "WWPS Soil % Clay 0-10cm", color="% Clay")
sam_map_clay

sam_map_sand <- ggmap(WWPS_map) + geom_point(data = sample_locations, aes(lon,lat,color=per_sand), size=2,alpha=0.7)+labs(x="longitude", y="latitude", title = "WWPS Soil % Sand 0-10cm", color="% Sand")
sam_map_sand

sam_map_silt <- ggmap(WWPS_map) + geom_point(data = sample_locations, aes(lon,lat,color=per_silt), size=2,alpha=0.7)+labs(x="longitude", y="latitude", title = "WWPS Soil % Silt 0-10cm", color="% Silt")
sam_map_silt

summary(sample_locations$ph1)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.670   5.683   5.965   5.966   6.207   6.960 

        summary(sample_locations$ec1)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    9.6    85.4   105.5   116.2   136.2   375.5 

        summary(sample_locations$gwc)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.07188 0.20471 0.26766 0.26942 0.34005 0.49295 

        summary(sample_locations$per_clay)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.523   6.923   9.866  10.348  13.168  32.303 

        summary(sample_locations$per_sand)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  12.12   16.53   19.81   20.57   24.34   33.62 

        summary(sample_locations$per_silt)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  41.57   64.73   70.08   69.08   74.57   84.75 

Now I can look at colinearity among all the samples I collected:

just_var <- sample_locations[c(5,6,10,20,27:28)]
pairs(just_var[1:6])

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(just_var[1:6],upper.panel=panel.cor ,lower.panel=panel.smooth)

Then, we use the data from the original file to create a SPDF class object so that data is clearly distinguished from the coordinates.

coordinates(sample_locations) <- ~ lon + lat

Error in `coordinates<-`(`*tmp*`, value = ~lon + lat) : 
  setting coordinates cannot be done on Spatial objects, where they have already been set

In order to perform kriging, we must first create a variogram model for each variable and then fitting a model to each variable’s variogram

lzn.vgm.ph <- variogram((ph1)~lon+lat, sample_locations)
lzn.fit.ph <- fit.variogram(lzn.vgm.ph, model=vgm(c("Exp","Sph", "Gau", "Mat")))

No convergence after 200 iterations: try different initial values?No convergence after 200 iterations: try different initial values?No convergence after 200 iterations: try different initial values?No convergence after 200 iterations: try different initial values?

plot(lzn.vgm.ph, lzn.fit.ph)

lzn.vgm.ec <- variogram((ec1)~lon+lat, sample_locations)
lzn.fit.ec <- fit.variogram(lzn.vgm.ec, model=vgm(c("Exp","Sph", "Gau", "Mat")))

singular model in variogram fit

[1] "a possible solution MIGHT be to scale semivariances and/or distances"

singular model in variogram fit

[1] "a possible solution MIGHT be to scale semivariances and/or distances"

singular model in variogram fit

[1] "a possible solution MIGHT be to scale semivariances and/or distances"

singular model in variogram fit

[1] "a possible solution MIGHT be to scale semivariances and/or distances"

plot(lzn.vgm.ec, lzn.fit.ec)

lzn.vgm.gwc <- variogram((gwc)~lon+lat, sample_locations, width=0.0002)
lzn.fit.gwc <- fit.variogram(lzn.vgm.gwc, model=vgm(c("Exp","Sph", "Gau", "Mat")))

singular model in variogram fit

plot(lzn.vgm.gwc, lzn.fit.gwc)

lzn.vgm.clay <- variogram((per_clay)~lon+lat, sample_locations, width=0.0001)
lzn.fit.clay <- fit.variogram(lzn.vgm.clay, model=vgm(c("Exp","Sph", "Gau", "Mat")))

singular model in variogram fit

plot(lzn.vgm.clay, lzn.fit.clay)

lzn.vgm.sand <- variogram((per_sand)~lon+lat, sample_locations)
lzn.fit.sand <- fit.variogram(lzn.vgm.sand, model=vgm(c("Exp","Sph", "Gau", "Mat")))
plot(lzn.vgm.sand, lzn.fit.sand)

lzn.vgm.silt <- variogram((per_silt)~lon+lat, sample_locations)
lzn.fit.silt <- fit.variogram(lzn.vgm.silt, model=vgm(c("Exp","Sph", "Gau", "Mat")))

No convergence after 200 iterations: try different initial values?

plot(lzn.vgm.silt, lzn.fit.silt, diff=TRUE)

lzn.fit.ph

  model     psill        range
1   Nug 0.0000000 0.000000e+00
2   Sph 0.1656692 8.652802e-05

lzn.fit.ec

  model     psill        range
1   Nug  714.1345 0.000000e+00
2   Sph 1602.3557 6.465045e-05

lzn.fit.gwc

  model       psill       range
1   Nug 0.002608152 0.000000000
2   Exp 0.003623269 0.001107374

lzn.fit.clay

  model     psill       range
1   Nug  4.129199 0.000000000
2   Sph 18.971155 0.000240231

lzn.fit.sand

  model     psill       range
1   Nug  7.888791 0.000000000
2   Gau 18.976200 0.000181561

lzn.fit.silt

  model    psill        range
1   Nug 10.73700 0.0000000000
2   Sph 55.78812 0.0003088612

Now that we fit the variogram models to the data, we will now use this information in kreging

bbox(sample_locations)

          min       max
lon -94.23546 -94.23402
lat  36.06563  36.06665

#Create a grid to estimate values over
x_range <- as.numeric(c(-94.23546, -94.23402))
y_range <- as.numeric(c(36.06563, 36.06665))
# create an empty grid of values ranging from the xmin-xmax, ymin-ymax
sample.grid <- expand.grid(x = seq(from = x_range[1],
                   to = x_range[2], 
                   length.out=30),
                   y = seq(from = y_range[1],                                           to = y_range[2], 
                       length.out=30))  # expand points to grid
class(sample.grid)

[1] "data.frame"

sample.grid$number=seq(from=1, to=900, length.out=900)
plot1 <- sample_locations %>% as.data.frame %>%
  ggplot(aes(lat, lon)) + geom_point(size=.25) + coord_equal() + 
ggtitle("Points with measurements")
# this is clearly gridded over the region of interest
plot2 <- sample.grid %>% as.data.frame %>%
  ggplot(aes(y, x)) + geom_point(size=.25) + coord_equal() + 
  ggtitle("Points at which to estimate")
library(gridExtra)
grid.arrange(plot1, plot2, ncol=2)

coordinates(sample.grid) <- ~ x + y

lzn.kriged.ph <- krige(ph1~1, sample_locations, sample.grid, model=lzn.fit.ph)

[using ordinary kriging]

lzn.kriged.ec <- krige(ec1~1, sample_locations, sample.grid, model=lzn.fit.ec)

[using ordinary kriging]

lzn.kriged.gwc <- krige(gwc~1, sample_locations, sample.grid, model=lzn.fit.gwc)

[using ordinary kriging]

lzn.kriged.clay <- krige(per_clay~1, sample_locations, sample.grid, model=lzn.fit.clay)

[using ordinary kriging]

lzn.kriged.sand <- krige(per_sand~1, sample_locations, sample.grid, model=lzn.fit.sand)

[using ordinary kriging]

lzn.kriged.silt <- krige(per_silt~1, sample_locations, sample.grid, model=lzn.fit.silt)

[using ordinary kriging]

a<- lzn.kriged.ph %>% as.data.frame %>%
  ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient(low = "yellow", high="red") +
  scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
  theme_bw() + labs(x="longitude", y="latitude", title = "WWPS Pred pH 0-10cm", fill="Pred pH")
b<- lzn.kriged.ec %>% as.data.frame %>%
  ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient(low = "yellow", high="red") +
  scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
  theme_bw() + labs(x="longitude", y="latitude", title = "WWPS Pred EC 0-10cm", fill="Pred EC (uS/cm")
c<- lzn.kriged.gwc %>% as.data.frame %>%
  ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient(low = "yellow", high="red") +
  scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
  theme_bw() + labs(x="longitude", y="latitude", title = "WWPS Pred GWC 0-10cm", fill="Pred GWC (% g/g)")
d<- lzn.kriged.clay %>% as.data.frame %>%
  ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient(low = "yellow", high="red") +
  scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
  theme_bw() + labs(x="longitude", y="latitude", title = "WWPS Pred % Clay 0-10cm", fill="Pred Clay (%)")
e<- lzn.kriged.sand %>% as.data.frame %>%
  ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient(low = "yellow", high="red") +
  scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
  theme_bw() + labs(x="longitude", y="latitude", title = "WWPS Pred % Sand 0-10cm", fill="Pred Sand (%)")
f<- lzn.kriged.silt %>% as.data.frame %>%
  ggplot(aes(x=x, y=y)) + geom_tile(aes(fill=var1.pred)) + coord_equal() +
  scale_fill_gradient(low = "yellow", high="red") +
  scale_x_continuous(labels=comma) + scale_y_continuous(labels=comma) +
  theme_bw() + labs(x="longitude", y="latitude", title = "WWPS Pred % Silt 0-10cm", fill="Pred Silt (%)")
a

b

c

d

e

f