library(sp)
library(rgdal)
## rgdal: version: 1.2-8, (SVN revision 663)
## 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/Jorn/Documents/R/win-library/3.4/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: C:/Users/Jorn/Documents/R/win-library/3.4/rgdal/proj
## Linking to sp version: 1.2-5
library(ggmap)
## Loading required package: ggplot2
library(ggplot2)
require(lattice)
## Loading required package: lattice
require(gstat)
## Loading required package: gstat
require(nlme)
## Loading required package: nlme
#metadata of the data set
#The data distributed to the SIC97 participants were 100 daily rainfall measurements made in Switzerland on the 8th of May 1986 which were randomly extracted from a dataset of 467 measurements. The participants had to estimate the rainfall at the 367 remaining locations. The measurements were in units of 1/10th of a mm, but values ranging from 1 to 10 were often observed due to air condensation. Figure 1 shows the locations of the 367 measurements for which the values had to be estimated (diamonds) while those of the 100 measurements used for the estimation are presented with proportional symbols (filled circles).
wd <- "C:/Users/Jorn/Downloads/spat-interp-comparison-1997"# setting the work space
setwd(wd)
swissRain = read.table(("sic_obs.dat"), sep=',',col.names=c('ID','x','y','rain'))
# Check the data set
print(swissRain)
## ID x y rain
## 1 13 -140463 -30977 151
## 2 14 -136211 -12166 255
## 3 22 -123184 -9111 79
## 4 23 -121276 9758 191
## 5 24 -120674 1960 194
## 6 29 -116780 -32596 334
## 7 30 -115480 -4816 107
## 8 35 -109680 20655 296
## 9 36 -109192 -36075 394
## 10 37 -109172 -34963 394
## 11 52 -93540 -14096 324
## 12 55 -92287 65956 105
## 13 66 -84732 68067 135
## 14 71 -83690 -18690 585
## 15 73 -83783 -80965 114
## 16 84 -75594 17901 334
## 17 95 -68707 -43345 131
## 18 102 -63854 -92327 78
## 19 105 -62677 19974 398
## 20 126 -46443 -37998 141
## 21 130 -45400 -1310 192
## 22 136 -42594 70948 151
## 23 138 -42513 82067 107
## 24 159 -34132 -26960 145
## 25 168 -29884 57526 334
## 26 172 -28425 47511 327
## 27 178 -26758 81974 213
## 28 181 -25311 69736 331
## 29 185 -24717 31927 400
## 30 188 -23978 27477 327
## 31 192 -21560 67497 380
## 32 198 -19610 -25914 94
## 33 202 -17218 8550 185
## 34 203 -17186 19669 239
## 35 207 -14028 71922 330
## 36 208 -13567 -57065 30
## 37 218 -6524 61902 254
## 38 224 254 -1483 53
## 39 226 1790 -51520 78
## 40 235 6322 16312 71
## 41 245 14022 -32601 62
## 42 246 13971 -10362 71
## 43 247 16178 18554 59
## 44 254 21414 38586 60
## 45 261 24364 56389 124
## 46 263 25738 86417 153
## 47 274 30369 60865 75
## 48 275 31263 34182 137
## 49 276 31216 43078 86
## 50 277 32301 -18075 129
## 51 281 33150 -33638 345
## 52 283 34145 -74775 441
## 53 287 33874 105361 184
## 54 292 37632 102049 121
## 55 293 39277 -34713 346
## 56 300 41380 -6899 270
## 57 302 42166 95408 100
## 58 305 44793 45389 45
## 59 314 45988 85429 107
## 60 317 47728 -39100 359
## 61 320 49386 -54654 278
## 62 322 49452 29858 72
## 63 335 54225 86611 131
## 64 336 56563 -83503 90
## 65 341 55862 72171 141
## 66 342 55851 73283 131
## 67 344 59158 -33438 452
## 68 357 63585 -91218 16
## 69 362 63172 91149 136
## 70 368 65718 64488 130
## 71 369 65456 87838 118
## 72 372 66256 83399 109
## 73 373 68438 24486 145
## 74 377 71226 -18851 254
## 75 378 70784 81227 140
## 76 381 71928 48992 152
## 77 384 73312 59018 60
## 78 392 78296 -34332 283
## 79 399 81859 40225 184
## 80 400 81960 86933 127
## 81 401 84694 230 220
## 82 406 85540 46949 178
## 83 408 87981 -16404 218
## 84 409 86885 58090 137
## 85 421 94901 74895 144
## 86 423 96720 9309 230
## 87 424 96467 24875 282
## 88 425 98187 -34040 129
## 89 436 105188 -40595 65
## 90 442 107070 -17206 190
## 91 449 113152 -15982 170
## 92 450 111636 64067 156
## 93 451 114546 -49321 131
## 94 455 120891 -21392 10
## 95 456 122870 -43597 99
## 96 458 125007 937 92
## 97 460 128096 -1222 67
## 98 466 142077 -45399 18
## 99 468 145302 -51996 20
## 100 471 150921 -61870 55
# the following seems to make the coordinates line up with epsg:2056
swissRain$x = swissRain$x - 17791.29 + 2672591
swissRain$y = swissRain$y - 13224.66 + 1200225
# the readme file says rain is in tenths of mm
swissRain$rain = swissRain$rain / 10
# Compare now the data set - how does it look now? - by printing
print(swissRain)
## ID x y rain
## 1 13 2514337 1156023 15.1
## 2 14 2518589 1174834 25.5
## 3 22 2531616 1177889 7.9
## 4 23 2533524 1196758 19.1
## 5 24 2534126 1188960 19.4
## 6 29 2538020 1154404 33.4
## 7 30 2539320 1182184 10.7
## 8 35 2545120 1207655 29.6
## 9 36 2545608 1150925 39.4
## 10 37 2545628 1152037 39.4
## 11 52 2561260 1172904 32.4
## 12 55 2562513 1252956 10.5
## 13 66 2570068 1255067 13.5
## 14 71 2571110 1168310 58.5
## 15 73 2571017 1106035 11.4
## 16 84 2579206 1204901 33.4
## 17 95 2586093 1143655 13.1
## 18 102 2590946 1094673 7.8
## 19 105 2592123 1206974 39.8
## 20 126 2608357 1149002 14.1
## 21 130 2609400 1185690 19.2
## 22 136 2612206 1257948 15.1
## 23 138 2612287 1269067 10.7
## 24 159 2620668 1160040 14.5
## 25 168 2624916 1244526 33.4
## 26 172 2626375 1234511 32.7
## 27 178 2628042 1268974 21.3
## 28 181 2629489 1256736 33.1
## 29 185 2630083 1218927 40.0
## 30 188 2630822 1214477 32.7
## 31 192 2633240 1254497 38.0
## 32 198 2635190 1161086 9.4
## 33 202 2637582 1195550 18.5
## 34 203 2637614 1206669 23.9
## 35 207 2640772 1258922 33.0
## 36 208 2641233 1129935 3.0
## 37 218 2648276 1248902 25.4
## 38 224 2655054 1185517 5.3
## 39 226 2656590 1135480 7.8
## 40 235 2661122 1203312 7.1
## 41 245 2668822 1154399 6.2
## 42 246 2668771 1176638 7.1
## 43 247 2670978 1205554 5.9
## 44 254 2676214 1225586 6.0
## 45 261 2679164 1243389 12.4
## 46 263 2680538 1273417 15.3
## 47 274 2685169 1247865 7.5
## 48 275 2686063 1221182 13.7
## 49 276 2686016 1230078 8.6
## 50 277 2687101 1168925 12.9
## 51 281 2687950 1153362 34.5
## 52 283 2688945 1112225 44.1
## 53 287 2688674 1292361 18.4
## 54 292 2692432 1289049 12.1
## 55 293 2694077 1152287 34.6
## 56 300 2696180 1180101 27.0
## 57 302 2696966 1282408 10.0
## 58 305 2699593 1232389 4.5
## 59 314 2700788 1272429 10.7
## 60 317 2702528 1147900 35.9
## 61 320 2704186 1132346 27.8
## 62 322 2704252 1216858 7.2
## 63 335 2709025 1273611 13.1
## 64 336 2711363 1103497 9.0
## 65 341 2710662 1259171 14.1
## 66 342 2710651 1260283 13.1
## 67 344 2713958 1153562 45.2
## 68 357 2718385 1095782 1.6
## 69 362 2717972 1278149 13.6
## 70 368 2720518 1251488 13.0
## 71 369 2720256 1274838 11.8
## 72 372 2721056 1270399 10.9
## 73 373 2723238 1211486 14.5
## 74 377 2726026 1168149 25.4
## 75 378 2725584 1268227 14.0
## 76 381 2726728 1235992 15.2
## 77 384 2728112 1246018 6.0
## 78 392 2733096 1152668 28.3
## 79 399 2736659 1227225 18.4
## 80 400 2736760 1273933 12.7
## 81 401 2739494 1187230 22.0
## 82 406 2740340 1233949 17.8
## 83 408 2742781 1170596 21.8
## 84 409 2741685 1245090 13.7
## 85 421 2749701 1261895 14.4
## 86 423 2751520 1196309 23.0
## 87 424 2751267 1211875 28.2
## 88 425 2752987 1152960 12.9
## 89 436 2759988 1146405 6.5
## 90 442 2761870 1169794 19.0
## 91 449 2767952 1171018 17.0
## 92 450 2766436 1251067 15.6
## 93 451 2769346 1137679 13.1
## 94 455 2775691 1165608 1.0
## 95 456 2777670 1143403 9.9
## 96 458 2779807 1187937 9.2
## 97 460 2782896 1185778 6.7
## 98 466 2796877 1141601 1.8
## 99 468 2800102 1135004 2.0
## 100 471 2805721 1125130 5.5
# Data set variables adn number of observations/rows
str(swissRain)
## 'data.frame': 100 obs. of 4 variables:
## $ ID : int 13 14 22 23 24 29 30 35 36 37 ...
## $ x : num 2514337 2518589 2531616 2533524 2534126 ...
## $ y : num 1156023 1174834 1177889 1196758 1188960 ...
## $ rain: num 15.1 25.5 7.9 19.1 19.4 33.4 10.7 29.6 39.4 39.4 ...
# Rainfall values statitsics
summary(swissRain$rain)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.000 9.775 14.100 18.015 25.425 58.500
# Rain Guage ID with with maximum observed rainfall
swissRain[which.max(swissRain$rain), ]
## ID x y rain
## 14 71 2571110 1168310 58.5
# Rain Guage ID with with minimum observed rainfall
swissRain[which.min(swissRain$rain), ]
## ID x y rain
## 94 455 2775691 1165608 1
# Quantile - frequency distribution
quantile(swissRain$rain)
## 0% 25% 50% 75% 100%
## 1.000 9.775 14.100 25.425 58.500
# frequency distribution
stem(swissRain$rain)
##
## The decimal point is 1 digit(s) to the right of the |
##
## 0 | 1222355666667777788889999
## 1 | 0011111122233333333444444445555556788899999
## 2 | 122345567888
## 3 | 023333333556899
## 4 | 0045
## 5 | 9
# Exploratory graphics
# Histogram (simple) - frequency distribution
hist(swissRain$rain)
# A more enhanced histogram
#The frequency histogram revealed that the rainfall data do not follow the normal distribution, even after applying a number of transformations.
h <- hist(swissRain$rain, breaks = seq(0, 60, by = 10), plot = F)
str(h)
## List of 6
## $ breaks : num [1:7] 0 10 20 30 40 50 60
## $ counts : int [1:6] 27 41 13 16 2 1
## $ density : num [1:6] 0.027 0.041 0.013 0.016 0.002 0.001
## $ mids : num [1:6] 5 15 25 35 45 55
## $ xname : chr "swissRain$rain"
## $ equidist: logi TRUE
## - attr(*, "class")= chr "histogram"
plot(h, col = cm.colors(length(h$mids))[length(h$count) - rank(h$count)
+ 1], ylim = c(0, max(h$count) + 10),
main = "Frequency histogram, Swiss Rainfall", sub = "Counts shown above bar, actual values shown with rug plot",
xlab = "Rainfall, mm")
rug(swissRain$rain) # adding the measurment on the botton with an small horizontal line
text(h$mids, h$count, h$count, pos = 3) # count adding on top of the bar
hlog10 <- hist((log10(swissRain$rain)), breaks = seq(0, 2, by = 0.1), plot = F)
str(hlog10)
## List of 6
## $ breaks : num [1:21] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
## $ counts : int [1:20] 1 0 2 1 1 0 1 6 9 6 ...
## $ density : num [1:20] 0.1 0 0.2 0.1 0.1 0 0.1 0.6 0.9 0.6 ...
## $ mids : num [1:20] 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 ...
## $ xname : chr "(log10(swissRain$rain))"
## $ equidist: logi TRUE
## - attr(*, "class")= chr "histogram"
plot(hlog10, col = heat.colors(length(hlog10$mids))[length(hlog10$count) - rank(hlog10$count) + 1], ylim = c(0, max(hlog10$count) + 5),
main = "Frequency histogram, Swiss Rainfall (log10)",
sub = "Counts shown above bar, actual values shown with rug plot",
xlab = "Log10 - Rainfall (mm)")
rug(log10(swissRain$rain))
text(hlog10$mids, hlog10$count, hlog10$count, pos = 3)
hsqrt <- hist((sqrt(swissRain$rain)), breaks = seq(0, 8, by = 0.5), plot = F)
str(hsqrt)
## List of 6
## $ breaks : num [1:17] 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 ...
## $ counts : int [1:16] 0 1 3 1 7 11 12 24 9 5 ...
## $ density : num [1:16] 0 0.02 0.06 0.02 0.14 0.22 0.24 0.48 0.18 0.1 ...
## $ mids : num [1:16] 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 ...
## $ xname : chr "(sqrt(swissRain$rain))"
## $ equidist: logi TRUE
## - attr(*, "class")= chr "histogram"
plot(hsqrt, col = topo.colors(length(hsqrt$mids))[length(hsqrt$count) - rank(hsqrt$count) + 1], ylim = c(0, max(hsqrt$count) + 5),
main = "Frequency histogram, Swiss Rainfall (Square Root)",
sub = "Counts shown above bar, actual values shown with rug plot",
xlab = "Square Root - Rainfall (mm)")
rug(sqrt(swissRain$rain))
text(hsqrt$mids, hsqrt$count, hsqrt$count, pos = 3)
hlog <- hist((log(swissRain$rain)), breaks = seq(0, 5, by = 0.5), plot = F)
str(hlog)
## List of 6
## $ breaks : num [1:11] 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 ...
## $ counts : int [1:10] 2 2 1 12 18 33 18 13 1 0
## $ density : num [1:10] 0.04 0.04 0.02 0.24 0.36 0.66 0.36 0.26 0.02 0
## $ mids : num [1:10] 0.25 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75
## $ xname : chr "(log(swissRain$rain))"
## $ equidist: logi TRUE
## - attr(*, "class")= chr "histogram"
plot(hlog, col = terrain.colors(length(hlog$mids))[length(hlog$count) - rank(hlog$count) + 1], ylim = c(0, max(hlog$count) + 5),
main = "Frequency histogram, Swiss Rainfall (log)",
sub = "Counts shown above bar, actual values shown with rug plot",
xlab = "Log - Rainfall (mm)")
rug(log(swissRain$rain))
text(hlog$mids, hlog$count, hlog$count, pos = 3)
# The frequency histogram revealed that the rainfall data do not follow the normal distribution. It reveal right-skewed distibution (this means high values are less). The transformation show a better distribution, especially the square root and log transformation.
# Visual spatial structure
coordinates(swissRain) <- c("x", "y") # MAKE A SPATIAL DATAFRAME BY DFINING X AND Y FIELDS
class(swissRain)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
print(bubble(swissRain, "rain", main = "Total Rainfall of the 8th of May 1986, Switzerland",
sub = "Rainfall (mm); symbol size proportional to value",
col = "limegreen", scales = list(draw = TRUE), maxsize = 2.0,
fill = F, xlab = "X (m)", ylab = "Y (m)", sizelab = "Y (m)",
xlim = c(2514000,2801000), ylim = c(1094000, 1293000), aspect = "iso",
panel = function(x,y, ...) {
panel.xyplot(x, y, ...)
panel.grid(h = -1, v = -1, col = "purple4", lty = 2)
},
key.entries = seq(0, 60, by = 10)))
# We tell R our data it's on LV95, EPSG=2056
swissRain@proj4string <- CRS("+init=epsg:2056")
# We can check the original coords
swissRain@coords
## x y
## [1,] 2514337 1156023
## [2,] 2518589 1174834
## [3,] 2531616 1177889
## [4,] 2533524 1196758
## [5,] 2534126 1188960
## [6,] 2538020 1154404
## [7,] 2539320 1182184
## [8,] 2545120 1207655
## [9,] 2545608 1150925
## [10,] 2545628 1152037
## [11,] 2561260 1172904
## [12,] 2562513 1252956
## [13,] 2570068 1255067
## [14,] 2571110 1168310
## [15,] 2571017 1106035
## [16,] 2579206 1204901
## [17,] 2586093 1143655
## [18,] 2590946 1094673
## [19,] 2592123 1206974
## [20,] 2608357 1149002
## [21,] 2609400 1185690
## [22,] 2612206 1257948
## [23,] 2612287 1269067
## [24,] 2620668 1160040
## [25,] 2624916 1244526
## [26,] 2626375 1234511
## [27,] 2628042 1268974
## [28,] 2629489 1256736
## [29,] 2630083 1218927
## [30,] 2630822 1214477
## [31,] 2633240 1254497
## [32,] 2635190 1161086
## [33,] 2637582 1195550
## [34,] 2637614 1206669
## [35,] 2640772 1258922
## [36,] 2641233 1129935
## [37,] 2648276 1248902
## [38,] 2655054 1185517
## [39,] 2656590 1135480
## [40,] 2661122 1203312
## [41,] 2668822 1154399
## [42,] 2668771 1176638
## [43,] 2670978 1205554
## [44,] 2676214 1225586
## [45,] 2679164 1243389
## [46,] 2680538 1273417
## [47,] 2685169 1247865
## [48,] 2686063 1221182
## [49,] 2686016 1230078
## [50,] 2687101 1168925
## [51,] 2687950 1153362
## [52,] 2688945 1112225
## [53,] 2688674 1292361
## [54,] 2692432 1289049
## [55,] 2694077 1152287
## [56,] 2696180 1180101
## [57,] 2696966 1282408
## [58,] 2699593 1232389
## [59,] 2700788 1272429
## [60,] 2702528 1147900
## [61,] 2704186 1132346
## [62,] 2704252 1216858
## [63,] 2709025 1273611
## [64,] 2711363 1103497
## [65,] 2710662 1259171
## [66,] 2710651 1260283
## [67,] 2713958 1153562
## [68,] 2718385 1095782
## [69,] 2717972 1278149
## [70,] 2720518 1251488
## [71,] 2720256 1274838
## [72,] 2721056 1270399
## [73,] 2723238 1211486
## [74,] 2726026 1168149
## [75,] 2725584 1268227
## [76,] 2726728 1235992
## [77,] 2728112 1246018
## [78,] 2733096 1152668
## [79,] 2736659 1227225
## [80,] 2736760 1273933
## [81,] 2739494 1187230
## [82,] 2740340 1233949
## [83,] 2742781 1170596
## [84,] 2741685 1245090
## [85,] 2749701 1261895
## [86,] 2751520 1196309
## [87,] 2751267 1211875
## [88,] 2752987 1152960
## [89,] 2759988 1146405
## [90,] 2761870 1169794
## [91,] 2767952 1171018
## [92,] 2766436 1251067
## [93,] 2769346 1137679
## [94,] 2775691 1165608
## [95,] 2777670 1143403
## [96,] 2779807 1187937
## [97,] 2782896 1185778
## [98,] 2796877 1141601
## [99,] 2800102 1135004
## [100,] 2805721 1125130
# We can plot the data to check
plot(swissRain)
# Then we transform the data to WGS84
swissRain2 <- spTransform(swissRain, CRS("+init=epsg:4326"))
# We can see the transformed coords
#swissRain2@coords
swissRain2
## coordinates ID rain
## 1 (6.321662, 46.55) 13 15.1
## 2 (6.373764, 46.71973) 14 25.5
## 3 (6.543686, 46.74867) 22 7.9
## 4 (6.565903, 46.91859) 23 19.1
## 5 (6.574934, 46.84851) 24 19.4
## 6 (6.63066, 46.53805) 29 33.4
## 7 (6.643931, 46.78805) 30 10.7
## 8 (6.716812, 47.01767) 35 29.6
## 9 (6.729983, 46.50742) 36 39.4
## 10 (6.730112, 46.51742) 37 39.4
## 11 (6.932047, 46.70622) 52 32.4
## 12 (6.941785, 47.42634) 55 10.5
## 13 (7.041771, 47.44571) 66 13.5
## 14 (7.061138, 46.66539) 71 58.5
## 15 (7.063823, 46.1052) 73 11.4
## 16 (7.165252, 46.99484) 84 33.4
## 17 (7.257662, 46.44409) 95 13.1
## 18 (7.321773, 46.00355) 102 7.8
## 19 (7.335033, 47.01377) 105 39.8
## 20 (7.547485, 46.49229) 126 14.1
## 21 (7.561815, 46.8223) 130 19.2
## 22 (7.60053, 47.4722) 136 15.1
## 23 (7.60191, 47.5722) 138 10.7
## 24 (7.708326, 46.59131) 159 14.5
## 25 (7.768379, 47.35112) 168 33.4
## 26 (7.787101, 47.26098) 172 32.7
## 27 (7.811284, 47.57088) 178 21.3
## 28 (7.829705, 47.46075) 181 33.1
## 29 (7.835053, 47.12065) 185 40.0
## 30 (7.844489, 47.08059) 188 32.7
## 31 (7.879284, 47.44043) 192 38.0
## 32 (7.897895, 46.60012) 198 9.4
## 33 (7.931937, 46.91) 202 18.5
## 34 (7.933276, 47.01001) 203 23.9
## 35 (7.979537, 47.47981) 207 33.0
## 36 (7.973988, 46.31956) 208 3.0
## 37 (8.078005, 47.38918) 218 25.4
## 38 (8.160056, 46.81853) 224 5.3
## 39 (8.17404, 46.36832) 226 7.8
## 40 (8.241952, 46.97807) 235 7.1
## 41 (8.335796, 46.53735) 245 6.2
## 42 (8.338455, 46.73739) 246 7.1
## 43 (8.371828, 46.99725) 247 5.9
## 44 (8.444036, 47.17684) 254 6.0
## 45 (8.486082, 47.33661) 261 12.4
## 46 (8.509683, 47.60649) 263 15.3
## 47 (8.566378, 47.37612) 274 7.5
## 48 (8.573107, 47.13603) 275 13.7
## 49 (8.574184, 47.21604) 276 8.6
## 50 (8.576797, 46.66588) 277 12.9
## 51 (8.584919, 46.52578) 281 34.5
## 52 (8.590018, 46.15564) 283 44.1
## 53 (8.621665, 47.77579) 287 18.4
## 54 (8.671101, 47.74549) 292 12.1
## 55 (8.664543, 46.51528) 293 34.6
## 56 (8.697758, 46.76515) 300 27.0
## 57 (8.730089, 47.68512) 302 10.0
## 58 (8.753904, 47.23492) 305 4.5
## 59 (8.778714, 47.59481) 314 10.7
## 60 (8.773674, 46.47458) 317 35.9
## 61 (8.791766, 46.33443) 320 27.8
## 62 (8.811832, 47.09452) 322 7.2
## 63 (8.888501, 47.60412) 335 13.1
## 64 (8.878084, 46.07381) 336 9.0
## 65 (8.906673, 47.474) 341 14.1
## 66 (8.906803, 47.484) 342 13.1
## 67 (8.923868, 46.52366) 344 45.2
## 68 (8.966883, 46.00322) 357 1.6
## 69 (9.008661, 47.64338) 362 13.6
## 70 (9.035314, 47.40318) 368 13.0
## 71 (9.038151, 47.61319) 369 11.8
## 72 (9.047578, 47.57313) 372 10.9
## 73 (9.060387, 47.04295) 373 14.5
## 74 (9.085098, 46.65268) 377 25.4
## 75 (9.107128, 47.55275) 378 14.0
## 76 (9.113182, 47.26267) 381 15.2
## 77 (9.134328, 47.35256) 384 6.0
## 78 (9.172964, 46.51209) 392 28.3
## 79 (9.241702, 47.18185) 399 18.4
## 80 (9.257318, 47.60182) 400 12.7
## 81 (9.266798, 46.8216) 401 22.0
## 82 (9.292347, 47.24155) 406 17.8
## 83 (9.304676, 46.67132) 408 21.8
## 84 (9.313623, 47.34144) 409 13.7
## 85 (9.42529, 47.49076) 421 14.4
## 86 (9.427359, 46.90061) 423 23.0
## 87 (9.429222, 47.04065) 424 28.2
## 88 (9.432149, 46.51046) 425 12.9
## 89 (9.521068, 46.44988) 436 6.5
## 90 (9.553759, 46.65975) 442 19.0
## 91 (9.633638, 46.66925) 449 17.0
## 92 (9.64323, 47.38939) 450 15.6
## 93 (9.639625, 46.36911) 451 13.1
## 94 (9.732658, 46.61861) 455 1.0
## 95 (9.749943, 46.41843) 456 9.9
## 96 (9.795125, 46.81829) 458 9.2
## 97 (9.83472, 46.79803) 460 6.7
## 98 (9.998873, 46.39686) 466 1.8
## 99 (10.03794, 46.3366) 468 2.0
## 100 (10.10651, 46.24613) 471 5.5
# Visual spatial structure
#coordinates(swissRain2) <- c("x", "y") # MAKE A SPATIAL DATAFRAME BY DFINING X AND Y FIELDS
class(swissRain2)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
# source of information: https://www.r-bloggers.com/google-maps-and-ggmap/
# example with roadmap background, 0.5 transparency, blue colour
mapImage <- get_map(location = c(lon = 8.5, lat = 46.5),
color = "color",
source = "google",
maptype = "roadmap",
zoom = 7)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=46.5,8.5&zoom=7&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
ggmap(mapImage) + geom_point(aes(x = x, y = y, size = swissRain2$rain), data = as.data.frame(coordinates(swissRain2)),
alpha = .5, color="blue")
# example with satellite background, 0.6 transparency, darkred colour
mapImage <- get_map(location = c(lon = 8.5, lat = 46.5),
color = "color",
source = "google",
maptype = "satellite",
zoom = 7)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=46.5,8.5&zoom=7&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false
ggmap(mapImage) + geom_point(aes(x = x, y = y, size = swissRain2$rain), data = as.data.frame(coordinates(swissRain2)),
alpha = .6, color="darkred")
# example with terrain background, 0.4 transparency, darkblue colour
mapImage <- get_map(location = c(lon = 8.5, lat = 46.5),
color = "color",
source = "google",
maptype = "terrain",
zoom = 7)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=46.5,8.5&zoom=7&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
ggmap(mapImage) + geom_point(aes(x = x, y = y, size = swissRain2$rain), data = as.data.frame(coordinates(swissRain2)),
alpha = .4, color="darkblue")
#Local spatial dependence
#So far the plot reveals a trend, but it does not show local spatial dependence. Therefore, a variogram cloud is generated, which shows the semivariances between all point pairs. This cloud reveals that semivariances generally increase with separation. However, this variogram shows too much details, as a result, an empirical variogram is created.
# Square root variogram cloud
vcsqrt <- variogram(sqrt(swissRain2$rain) ~ 1, loc = swissRain2, cloud = T)
plot(vcsqrt, main = "Variogram cloud, SQRT(Rainfall)", xlab = "separation distance, km")
# Log variogram cloud
vclog <- variogram(log(swissRain2$rain) ~ 1, loc = swissRain2, cloud = T)
plot(vclog, main = "Variogram cloud, LOG(Rainfall)", xlab = "separation distance, km")
#The empirical variogram
# The variogram cloud shows all point-pairs, but is clearly far too detailed to be useful for examining the general pattern of spatial dependence. For that we use the empirical variogram, which organizes the cloud into bins, like a histogram
# Square root empirical variogram
vsqrt <- variogram(sqrt(swissRain2$rain) ~ 1, loc = swissRain2)
print(vsqrt)
## np dist gamma dir.hor dir.ver id
## 1 14 4.882604 0.08444402 0 0 var1
## 2 64 11.574312 0.50269273 0 0 var1
## 3 109 19.290215 0.56620519 0 0 var1
## 4 137 27.468503 1.13546978 0 0 var1
## 5 135 35.233659 1.17946862 0 0 var1
## 6 187 42.527917 1.41172138 0 0 var1
## 7 175 50.424014 1.85932231 0 0 var1
## 8 207 58.163794 2.09414479 0 0 var1
## 9 225 65.801721 2.01550982 0 0 var1
## 10 227 73.878875 2.21856712 0 0 var1
## 11 228 81.484852 2.22385899 0 0 var1
## 12 230 89.482686 2.04107113 0 0 var1
## 13 245 96.842408 2.20268283 0 0 var1
## 14 276 104.963383 1.82401282 0 0 var1
## 15 266 112.509861 1.56930786 0 0 var1
print(plot(vsqrt, pl = T, main = "Empirical variogram, SQRT(Rain)",
xlab = "separation distance, km", col = "darkblue",
pch = 20))
# Log empirical variogram
vlog <- variogram(log(swissRain2$rain) ~ 1, loc = swissRain2)
print(vlog)
## np dist gamma dir.hor dir.ver id
## 1 14 4.882604 0.02527802 0 0 var1
## 2 64 11.574312 0.24611211 0 0 var1
## 3 109 19.290215 0.21767310 0 0 var1
## 4 137 27.468503 0.34368105 0 0 var1
## 5 135 35.233659 0.35049541 0 0 var1
## 6 187 42.527917 0.43663970 0 0 var1
## 7 175 50.424014 0.52989186 0 0 var1
## 8 207 58.163794 0.57109265 0 0 var1
## 9 225 65.801721 0.58066497 0 0 var1
## 10 227 73.878875 0.63704408 0 0 var1
## 11 228 81.484852 0.64700329 0 0 var1
## 12 230 89.482686 0.61750686 0 0 var1
## 13 245 96.842408 0.58297835 0 0 var1
## 14 276 104.963383 0.52919536 0 0 var1
## 15 266 112.509861 0.51044998 0 0 var1
print(plot(vlog, pl = T, main = "Empirical variogram, Log(Rain)",
xlab = "separation distance, km", col = "darkblue",
pch = 20))
# Visualizing anisotropy
# The auto-correlated spatial process that produced the dependence may be stronger in one direction than in others. This type of spatial dependence is called anisotropic (Rossiter & Truong, 2016b). Two visualization methods are applied to detect anisotropy. Firstly, a variogram map is generated. Secondly, directional variograms is created.
# Variogram Maps
print(plot(variogram(sqrt(rain) ~ 1, swissRain2, map = TRUE,
cutoff = 80, width = 10), main = "Variogram map, Rainfall Switzerland, SQRT Rainfall",
sub = "Cutoff of 80 km and bin width of 10 km",
col.regions = bpy.colors(64)))
#anisotropic: 180° (sOUTH)
print(plot(variogram(log(rain) ~ 1, swissRain2, map = TRUE,
cutoff = 80, width = 10), main = "Variogram map, Rainfall Switzerland, LOG Rainfall",
sub = "Cutoff of 80 km and bin width of 10 km",
col.regions = bpy.colors(64)))
# anisotropic: 180° (SOUTH)
# However, the autocorrelated spatial process that produced the dependence may be stronger in one direction than in others. This type of spatial dependence is called anisotropic, from the Greek "an-" (not) + isotropic. In practice, we consider the case where there is one direction of strongest spatial dependence, with the weakest spatial dependence at right angles (orthogonal) to it.
# Directional Variograms