The properties of the used data can affect the spatial estimation. Therefore we determine the statistical and spatial statistical characteristics in this first part of the seminar. I refer for the whole of process of geostatistics to a course run by Faculty of Geo-Information Science and Earth Observation (University of Twente, https://www.itc.nl/C18-GFM-DE-05) (Rossiter & Truong).

Content

1. Data Prep
1.1 Metadata of the Data Set
1.2 Data Set and some Prep
2. Data visualisation and Exploration
2.1 Descriptive Statistics
2.2 Exploratory Graphics
2.3 Visual Spatial Structure
- 2.3.1 Basic Postplots
- 2.3.2 Postplots with Topographical Features
2.4 Local spatial dependence
- 2.4.1 Variogram Cloud
- 2.4.2 Empirical Variogram Variogram Cloud
2.5 Visualizing Anisotropy
- 2.5.1 Variogram Maps
- 2.5.2 Directional Variograms
References

Required R modules

require(sp)

## Loading required package: sp

require(rgdal)

## Loading required package: 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

require(ggmap)

## Loading required package: ggmap

## Loading required package: ggplot2

require(ggplot2)
require(lattice)

## Loading required package: lattice

require(gstat)

## Loading required package: gstat

require(nlme)

## Loading required package: nlme

1. Data Prep

1.1 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.

Source Data Set: http://mldata.org/repository/data/viewslug/spat-interp-comparison-1997/1/

1.2 Data Set and some Prep

Assigning some field names to the data set.

swissRain = read.table(("C:/Users/Jorn/Downloads/spat-interp-comparison-1997/sic_obs.dat"), sep=',',col.names=c('ID','x','y','rain'))

Check the data set by a print statement

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

Lining up the coordinates with epsg:2056 (http://spatialreference.org/ref/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 - converting to 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

2.Data visualisation and Exploration

2.1 Descriptive Statistics

Print the data set variables and 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

2.2 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.

Normal Histogram

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

Log10 Histogram

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)

Square Root Histogram

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)

Log Histogram

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 transformations show a better distribution, especially the square root and log transformation.

2.3 Visual Spatial Structure

Making a spatial dataframe by defining x and y fields

coordinates(swissRain) <- c("x", "y") 
class(swissRain)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

2.3.1 Basic Postplots

The most basic visualization of spatial distribution in two dimensions is the postplot. It shows the location of each observation point in geographic space, with the value of the attribute in feature space shown by the symbol. For continuous variables, this is generally by making the symbol size proportional to the attribute value.

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)))

2.3.2 Postplots with Topographical Features

Telling R our data it’s on LV95, EPSG=2056

swissRain@proj4string <- CRS("+init=epsg:2056")

Checking 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

Then we transform the data to WGS84

swissRain2 <- spTransform(swissRain, CRS("+init=epsg:4326"))

We can see the transformed 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

Making a spatial dataframe by defining x and y coordinates

class(swissRain2)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

source of information of the backgrounds : https://www.r-bloggers.com/google-maps-and-ggmap/ The first 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")

2.4 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.

2.4.1 Variogram Cloud

The 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")

2.4.2 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

Some important Variogram terms:

nugget: the value where the semivariogram model intercepts the y-axis.
range: The distance where the model first flattens out is known as the range.
sill: The value that the semivariogram model attains at the range (the value on the y-axis) is called the sill.
partial sill: The partial sill is the sill minus the nugget.

The 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))

This variogram:

nugget: <0
range: 70
sill: 2.2

The 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))

This variogram:

nugget:<0
range: 75
sill: 0.6

2.5 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). Two visualization methods are applied to detect anisotropy. Firstly, a variogram map is generated. Secondly, directional variograms is created.

2.5.1 Variogram Maps

There are several pre-defined methods that create colour ramps:

heat.colors: from red to white-hot
terrain.colors
topo.colors: as used in topographic maps
cm.colors: cyan through magenta (pastels)
bpy.colors: blue-purple-yellow (looks good printed)
hsv: custom ramp specifying hue, saturation and value
rainbow: custom ramp from a given start to stop colour
gray: gray scale

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: 150 - 210°

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: 50 - 210°

The cutoff specified with the cutoff optional argument, must be an integer multiple of the bin width, specified with the width optional argument, in order to obtain a (correct) symmetric variogram map. In the above example cutoff=80, which is exactly eight times the width=10.

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.

2.5.2 Directional Variograms