setwd("/Users/shannoncall/Desktop/ESCI597A/Week06")
load("wUSPrec_Bnd.Rdata")
require(gstat)
## Loading required package: gstat
require(spatstat)
## Loading required package: spatstat
## Loading required package: nlme
##
## spatstat 1.45-0 (nickname: 'One Lakh')
## For an introduction to spatstat, type 'beginner'
##
## Note: spatstat version 1.45-0 is out of date by more than 12 weeks; a newer version should be available.
##
## Attaching package: 'spatstat'
## The following object is masked from 'package:gstat':
##
## idw
require(sp)
## Loading required package: sp
require(ncf)
## Loading required package: ncf
require(spdep)
## Loading required package: spdep
## Loading required package: Matrix
require(raster)
## Loading required package: raster
##
## Attaching package: 'raster'
## The following objects are masked from 'package:spatstat':
##
## area, rotate, shift
## The following object is masked from 'package:nlme':
##
## getData
require(ggplot2)
## Loading required package: ggplot2
sp_map <- fortify(wUSBoundary)
## Regions defined for each Polygons
wUSPrec.df <- data.frame(wUSPrec)
wUSPrec.df
## prec long lat optional
## 1 98.50996 -2.728566e+05 528678.350 TRUE
## 2 84.43426 -2.512029e+05 516639.677 TRUE
## 3 94.24108 -2.654868e+05 528352.447 TRUE
## 4 81.44355 -2.507462e+05 527727.412 TRUE
## 63 80.24970 4.401966e+04 556007.247 TRUE
## 68 128.07258 1.100441e+05 556846.660 TRUE
## 840 26.70616 -4.842403e+04 -401100.314 TRUE
## 1049 40.85034 3.549610e+05 -342265.575 TRUE
## 1051 31.82195 1.853417e+05 -341072.916 TRUE
## 1053 26.32022 1.303850e+05 -422252.612 TRUE
## 1054 37.07461 2.631869e+05 -370783.307 TRUE
## 1055 32.44568 2.112818e+05 -334778.510 TRUE
## 1056 34.91593 -1.417950e+04 -273547.329 TRUE
## 1061 30.99387 -1.950007e+05 -385233.780 TRUE
## 1062 30.27404 -2.483327e+05 -357996.815 TRUE
## 1063 45.82434 -1.941697e+05 -295276.858 TRUE
## 1066 34.62001 -2.963196e+05 -293874.642 TRUE
## 1067 47.85421 -2.195228e+05 -250051.430 TRUE
## 1068 35.06312 -2.798315e+05 -238959.759 TRUE
## 1070 40.87743 -1.476765e+05 -197495.565 TRUE
## 1071 46.39114 -2.243427e+05 -194289.730 TRUE
## 1072 34.02764 -2.752868e+05 -167940.581 TRUE
## 1073 64.35465 -9.041479e+04 -172922.134 TRUE
## 1074 60.92280 -2.060883e+05 -137019.936 TRUE
## 1075 99.13415 -1.745293e+05 -118996.395 TRUE
## 1076 40.18229 -2.729419e+05 -110160.393 TRUE
## 1077 58.15568 -9.709868e+04 -67160.331 TRUE
## 1195 45.11766 2.168176e+05 -114429.193 TRUE
## 1196 41.67888 3.288043e+05 -35268.020 TRUE
## 1197 35.43456 2.665202e+05 -31374.474 TRUE
## 1199 37.90724 2.202623e+05 2554.619 TRUE
## 1200 30.65565 2.446690e+05 14546.780 TRUE
## 1202 55.15561 2.625081e+05 68683.952 TRUE
## 1203 85.38318 2.922124e+05 114495.486 TRUE
## 1204 47.59461 -3.052262e+04 -197915.802 TRUE
## 1205 49.44257 -4.419014e+04 -143355.994 TRUE
## 1206 44.56676 9.401851e+04 -78329.700 TRUE
## 1207 65.55674 1.721239e+05 95737.991 TRUE
## 1208 74.15409 1.655181e+05 137862.689 TRUE
## 1209 55.38613 -1.490096e+04 137990.769 TRUE
## 1210 50.79841 3.521281e+04 152533.832 TRUE
## 1211 55.88187 9.214237e+04 166447.389 TRUE
## 1212 29.18475 5.668004e+04 202753.182 TRUE
## 1213 44.78784 -2.733908e+05 -81198.327 TRUE
## 1215 60.34926 -2.673250e+05 -31342.156 TRUE
## 1216 70.74156 -2.467086e+05 -18779.268 TRUE
## 1217 45.53637 -1.031989e+05 8594.478 TRUE
## 1218 117.11843 -1.692438e+05 23312.746 TRUE
## 1219 126.54753 -1.980355e+05 48605.160 TRUE
## 1221 49.57582 -7.244522e+04 40455.348 TRUE
## 1222 62.09862 -2.534870e+05 76123.283 TRUE
## 1223 54.44727 -2.486354e+05 140520.363 TRUE
## 1224 124.79335 -1.626598e+05 126660.134 TRUE
## 1225 182.86458 -1.687412e+05 164658.007 TRUE
## 1226 32.66142 -8.906517e+04 163066.063 TRUE
## 1227 65.68795 -2.418116e+05 175890.622 TRUE
## 1228 85.40244 -2.069654e+05 191373.340 TRUE
## 1229 37.66177 -1.182797e+05 189128.677 TRUE
## 1230 109.41849 -2.253390e+05 244307.710 TRUE
## 1300 102.75003 3.510343e+05 386968.275 TRUE
## 1310 97.24472 3.305902e+05 499307.801 TRUE
## 1313 31.20381 -5.622367e+04 166034.226 TRUE
## 1314 35.21251 -6.366777e+04 202818.753 TRUE
## 1315 31.95533 -4.185270e+04 330538.567 TRUE
## 1316 47.18121 -7.200867e+04 355269.353 TRUE
## 1317 32.59864 -2.415699e+04 380464.743 TRUE
## 1318 65.36126 2.290445e+05 270037.906 TRUE
## 1319 148.59675 3.436154e+05 244021.410 TRUE
## 1320 65.89162 1.849234e+05 279724.935 TRUE
## 1321 87.35635 2.837530e+05 285607.833 TRUE
## 1322 95.97598 2.314582e+05 309073.249 TRUE
## 1323 70.44470 2.145460e+05 311812.893 TRUE
## 1324 72.19817 1.988595e+05 324647.395 TRUE
## 1325 24.32664 6.939418e+04 247375.768 TRUE
## 1326 23.68758 1.958497e-09 258120.087 TRUE
## 1327 125.97700 2.921703e+05 400616.850 TRUE
## 1328 59.59675 1.849518e+05 407618.681 TRUE
## 1329 110.24043 2.385044e+05 423899.064 TRUE
## 1330 122.64080 2.548467e+05 482350.312 TRUE
## 1331 40.07188 1.237000e+05 348294.914 TRUE
## 1332 35.30602 9.830887e+04 370078.066 TRUE
## 1333 52.37343 1.398004e+05 408660.452 TRUE
## 1334 55.66220 9.979007e+04 418999.583 TRUE
## 1335 120.78646 2.344988e+05 489315.666 TRUE
## 1336 116.05539 1.625993e+05 512524.341 TRUE
## 1337 117.86482 1.623586e+05 549145.983 TRUE
## 1339 78.27514 1.918065e+05 177546.424 TRUE
## 1340 53.07813 1.223735e+05 219247.367 TRUE
## 1341 71.33377 1.320954e+05 236129.866 TRUE
## 1342 58.51237 1.539280e+05 259983.879 TRUE
## 1343 47.89019 -5.261751e+03 407088.450 TRUE
## 1344 56.39172 -5.344248e+04 485062.092 TRUE
## 1345 57.74608 -1.407495e+04 497056.078 TRUE
## 1346 63.95402 1.849278e+04 505947.295 TRUE
## 1347 124.80247 -3.001463e+05 247393.497 TRUE
## 1348 161.21033 -3.022518e+05 169533.905 TRUE
## 1350 154.48627 -2.845347e+05 302340.359 TRUE
## 1351 96.50162 -2.254292e+05 306629.574 TRUE
## 1352 149.01293 -2.905845e+05 337120.337 TRUE
## 1353 272.76229 -1.313442e+05 381801.528 TRUE
## 1354 154.25992 -1.393227e+05 396424.457 TRUE
## 1355 97.01533 -1.728638e+05 409493.289 TRUE
## 1356 119.21800 -1.636265e+05 445913.737 TRUE
## 1368 74.63322 -2.049725e+05 461594.292 TRUE
## 1369 150.15614 -1.657957e+05 503728.132 TRUE
## 1370 84.74396 -2.068165e+05 517194.016 TRUE
## 1371 108.38844 -1.858031e+05 527619.585 TRUE
## 1372 152.65415 -1.717495e+05 554953.001 TRUE
## 1374 185.55245 -1.389987e+05 308582.494 TRUE
## 1375 150.70730 -1.516396e+05 354468.211 TRUE
## 1376 97.16667 -1.772859e+05 359589.854 TRUE
## 1377 97.62209 -2.170782e+05 370853.516 TRUE
## 1403 29.03150 3.813304e+05 -107114.263 TRUE
Bubble Plot of Precipitation
bubble(wUSPrec, xcol="long",ycol="lat",zcol="prec",maxsize=2.5,main="")
Plot of Precipitation over polygon
gg <- ggplot()
gg <- gg + geom_map(data=sp_map,map=sp_map,
aes(x=long, y=lat, map_id=id),
fill= "white",color="black")
gg <- gg + geom_point(data=wUSPrec.df, aes(x=long, y=lat, size=prec),
color="magenta")
gg <- gg + scale_size(breaks=seq(0,300,by=50))
gg <- gg + labs(x=NULL, y=NULL)
gg <- gg + theme(panel.grid=element_blank(), panel.border=element_blank())
gg <- gg + theme(axis.ticks=element_blank(), axis.text=element_blank())
gg <- gg + theme(legend.position="top")
gg
Scatter Plots
hscat(prec~1,data=wUSPrec,breaks = seq(0,1e5,by=2e4))
Precipitation is autocorrelated at near distances.
Moran’s I
w <- 1/as.matrix(dist(coordinates(wUSPrec)))
diag(w) <- 0
moran.test(wUSPrec$prec,mat2listw(w))
##
## Moran I test under randomisation
##
## data: wUSPrec$prec
## weights: mat2listw(w)
##
## Moran I statistic standard deviate = 11.886, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.1657032288 -0.0089285714 0.0002158469
Moran’s I is 0.17, indicating autocorrelation. We could adjust the default to include more near neighbors.
Moran’s I with Eight Nearest Neighbor
w8 <- knn2nb(knearneigh(wUSPrec,k=8))
moran.test(wUSPrec$prec,nb2listw(w8))
##
## Moran I test under randomisation
##
## data: wUSPrec$prec
## weights: nb2listw(w8)
##
## Moran I statistic standard deviate = 10.687, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.444845197 -0.008928571 0.001802824
YIKES! It’s worse! 0.44…. lets try 16 neighbors…
w16 <- knn2nb(knearneigh(wUSPrec,k=16))
moran.test(wUSPrec$prec,nb2listw(w16))
##
## Moran I test under randomisation
##
## data: wUSPrec$prec
## weights: nb2listw(w16)
##
## Moran I statistic standard deviate = 13.173, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.3675004778 -0.0089285714 0.0008166128
n <- 6
res <- data.frame(k=2^(1:n),I=rep(NA,n))
for(i in 1:n){
w <- knn2nb(knearneigh(wUSPrec,k=2^i))
res$I[i] <- moran.test(wUSPrec$prec,nb2listw(w))$estimate[1]
}
plot(res,type="b",main="Moran's I by number of neighbors",pch=20,cex=1.5)
As k increases, I decreases. This means the more neighbors we use, the lower I value we get. This doesn’t help out a bunch though, since the residuals are most likely not clean.
anivar=variogram(prec~1, data=wUSPrec,cloud=FALSE, alpha=c(0,90,180,270))
plot(anivar,pch=20,cex=1.25,main="Precipitation by direction")
This makes sense that weather is direction, especially since North to South direction follows along the Cascade mountain range. There is also an east/west direction, again, indicating a mountain range which influences weather patterns.
Variogram of Precipitaton
PrecVarCloud <- variogram(prec~1, wUSPrec, cloud = TRUE)
plot(PrecVarCloud,pch=20,cex=1.5,col="dodgerblue4",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)",main = "Precipitation")
Variogram of Precipitaton without the cloud
PrecVar <- variogram(prec~1, wUSPrec, cloud = FALSE)
plot(PrecVar,pch=20,cex=1.5,col="goldenrod3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (m)", main = "Precipitation")
High autocorrelation at distanced around 2e+05 meters… this makes sense. Weather is similar to near weather…
Spline
spline.Prec<-spline.correlog(x=coordinates(wUSPrec)[,1],
y=coordinates(wUSPrec)[,2],z=wUSPrec$prec, resamp=100,
quiet=TRUE)
plot(spline.Prec)
Precipitation is similar to precipiation near itself. This makes sense since the eastern side of Washington is much more dry than western Washington. We’d expect precipitation to be similar and autocorrelated with weather near it.