library(sp)
## Warning: package 'sp' was built under R version 3.3.3
library(maptools)
## Warning: package 'maptools' was built under R version 3.3.3
## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.3.3
## rgdal: version: 1.2-7, (SVN revision 660)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/Nattybug/Documents/R/win-library/3.3/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/Nattybug/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-4
library(raster)
## Warning: package 'raster' was built under R version 3.3.3
library(gstat)
## Warning: package 'gstat' was built under R version 3.3.3
data(meuse)
head(meuse)
##        x      y cadmium copper lead zinc  elev       dist   om ffreq soil
## 1 181072 333611    11.7     85  299 1022 7.909 0.00135803 13.6     1    1
## 2 181025 333558     8.6     81  277 1141 6.983 0.01222430 14.0     1    1
## 3 181165 333537     6.5     68  199  640 7.800 0.10302900 13.0     1    1
## 4 181298 333484     2.6     81  116  257 7.655 0.19009400  8.0     1    2
## 5 181307 333330     2.8     48  117  269 7.480 0.27709000  8.7     1    2
## 6 181390 333260     3.0     61  137  281 7.791 0.36406700  7.8     1    2
##   lime landuse dist.m
## 1    1      Ah     50
## 2    1      Ah     30
## 3    1      Ah    150
## 4    0      Ga    270
## 5    0      Ah    380
## 6    0      Ga    470
class(meuse)
## [1] "data.frame"
coordinates(meuse)<-~x+y

proj4string(meuse)<-CRS("+init=epsg:28992")
class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
bubble(meuse,xcol="x",ycol="y",zcol="cadmium",maxsize=2.5,main="Cadmium concentrations (ppm)")

##can use the get data function in the raster package to find more data sets
cadvarCloud<-variogram(cadmium~1,meuse,cloud=TRUE)
plot(cadvarCloud,pch=20,cex=1.5,col="black",ylab=expression("Semivariance("*gamma*")"),xlab="Distance (m)",main="Cadmium concentrations (ppm)")

## sp package use data(meuse.rv) to get river outline map
##rasterr get data can give you borderlines, elevation, or climate
##mapmisc gives you GNcities
library(mapmisc)
## Warning: package 'mapmisc' was built under R version 3.3.3
try(CRS("init=espg:3857"),silent=TRUE)
library(sp)
library(maptools)
map.new(meuse)

library(sp)
library(raster)
map.new(meuse)

library(mapmisc)
map.new(meuse, legendRight=FALSE,buffer=0,mar=c(0,0,0,0))
if(require('rgdal',quietly=TRUE)) {proj4string(meuse)<-CRS("+init=epsg:28992")} else{proj4string(meuse)<-CRS(paste("+proj=sterea+lat_0=52.1561605555555+lon_0+5.3876888888888","+k-0.9999079+x_0=155000+y_0=463000+ellps=bessel+units=m+no_defs"))}
meuse$soilFac=factor(meuse$soil,levels=c(1,2,3),labels=c("calcareous","non-calcareous","redbrick"))

##I'm going to try looking at soil type for the meuse set
soilcol=colourScale(meuse$soilFac)
data("netherlands")
map.new(meuse)
plot(nldTiles,add=TRUE)
points(nldCities)
text(nldCities,label=nldCities$name,pos=2)
points(meuse,pch=16,col=soilcol$plot)
legend('topleft',fill=soilcol$col,legend=soilcol$legend)
insetMap(meuse,"bottomright",map=world)

Next I’m going to use the copper data as my z factor and look at a semivariogram versus a log transformed z variable variogram

copvarCloud<-variogram(copper~1,meuse,cloud=TRUE)
copvarCloud2<-variogram(copper~1,meuse,cloud=FALSE)
plot(copvarCloud,pch=20,cex=1.5,col="black",ylab=expression("Semivariance("*gamma*")"),xlab="Distance (m)",main="Copper concentrations (ppm)")

plot(copvarCloud2,pch=20,cex=1.5,col="black",ylab=expression("Semivariance("*gamma*")"),xlab="Distance (m)",main="Copper concentrations (ppm)")

copvarCloudlog<-variogram(log(copper)~1,meuse,cloud=TRUE)
plot(variogram(log(copper)~1, meuse),main="Log Transformed Copper Concentrations",ylab=expression("Semivariance("*gamma*")"),xlab="Distance in m")

It appears that the sill shows up at the same place in both graphs and the overall shape looks very similar so perhaps copper is more normally distributed than I thought.The sill appears at about 600m, so once the points are 600m apart or farther they are least similar. The y acis was different of course becaause of the log transformation. The map and width arguments are tied to how much of the map shows and the width of the cell size.

plot(variogram(log(copper)~1, meuse,cutoff=1300,width=100))

plot(variogram(log(copper)~1, meuse,cutoff=500))

plot(variogram(log(copper)~1, meuse,cutoff=1300,width=50))

If the width is set to a smaller number, the variogram because more “noisy” and less smooth. Using a larger width is almost like applying a filter. Changing the cutoff seems dangerous because if the cutoff is set too low, you may miss the sill!

load("~/Western/ESCI 597/mysteryData.Rdata")
coordinates(dat)=c("x","y")
plot(variogram(z~1, dat),xlab="distance in angstroms",main="Keebler elf concentration")

the sill is hard to see so I will change the cutoff argument and see if it becomes clear

plot(variogram(z~1, dat,cutoff=30),xlab="distance in angstroms",main="keebler elf concentration")

Oh Lordy what the heck is going on here. It seem as though the sill was at 15 angstroms but then starts trending down again. This makes me think that if the points are less than a distance of 15 and more than a distance of 20 apart, they are more similar to eachother than if they are 15-20 meters apart.It may not be appropriate to compare points that far out so I will say the sill is at a distance of 15 and stop there.

plot(variogram(z~1, dat, width=5),xlab="distance in angstroms",main="Keebler elf concentration")

Changing the width to a larger bin and using the default cutoff makes the sill at 15 more pronounced.