-
Notifications
You must be signed in to change notification settings - Fork 2
/
index.html
1109 lines (901 loc) · 50.4 KB
/
index.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="" xml:lang="">
<head>
<title>Working with spatial data in R</title>
<meta charset="utf-8" />
<meta name="date" content="2019-05-22" />
<link href="libs/remark-css/default.css" rel="stylesheet" />
<link href="libs/tabwid/tabwid.css" rel="stylesheet" />
<script src="libs/tabwid/tabwid.js"></script>
<link rel="stylesheet" href="custom-fonts.css" type="text/css" />
</head>
<body>
<textarea id="source">
class: center, middle, inverse, title-slide
# Working with spatial data in R
## <a href="https://github.com/ecodatasci-tlv/spatial" class="uri">https://github.com/ecodatasci-tlv/spatial</a>
### 2019-05-22
---
# Topics we'll cover today:
* Types of spatial data
* Importing/Exporting spatial data
* Projections, extents, and units
* Comparison of different GIS tools out there
* Spatial data manipulation
* Intro to spatial analysis
---
class: inverse, middle, center
#Data types
.small[
https://www.gislounge.com/geodatabases-explored-vector-and-raster-data/
]
---
.center[
.img-med[
]]
---
class: middle, center
# **Vector data**
Equivalent to a spreadsheet with a geometry column that describes how to plot the data that represents that feature
---
### **Points**
Used to represent nonadjacent features
.center[
.img-small[
]]
---
### **Lines**
One dimensional and are only used to measure length
.center[
.img-small[
]]
---
### **Polygons**
Two dimensional and are used for measurement of area and parameters of geographical feature.
.center[
.img-small[
]]
---
# File type
* Polygon and lines come in a shapefile format, point usually are stored in XY table (a spreadsheet with latitude and longitude columns)
* A shapefile is structured from 4 main files:
- **shp** - the geometry of each feature
- **dbf** - a dBase file with the attribute data for each feature (can be opened in excel; stores the metadata of each feature)
- **prj** - The projection file
- **shx** - a spatial index that helps a GIS program find the features in the .shp file quickly
**IMPORTANT:** each shapefile can only store one type of geometry - point/ polygon/ line
---
class: middle, center
# **Raster data**
Cell-based data, represents surfaces. Each cell represents a value
---
## **Raster can represent two types of data**
### **Continuous** - temperature, elevation, etc.
### **Discrete** - population density, species richness, etc.
---
# **Rasters have three types of datasets:**
### **Thematic**
Used to represent discrete data, such as soil type or land use
.center[
.img-small[
]]
---
### **Spectral**
Aerial or setalite imagery used during digital map generation
.center[
.img-med[
]]
---
### **Surface**
Represents continuous changes across a landscape (for example elevation)
.center[
.img-med[
]]
---
class: center
# Vector and raster comparison
<div class="tabwid"><table style='border-collapse:collapse;'><thead><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:bold;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Data_type</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:bold;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Pros</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:bold;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Cons</span></p></td></tr></thead><tbody><tr><td rowspan="3" style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:bold;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Vector</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">High geographical accuracy</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Doesn’t work well with continuous data</span></p></td></tr><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Follows topology rules which increases integrity</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Processing intensive due to topology checks</span></p></td></tr><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Used for measurements of proximity</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;"></span></p></td></tr><tr><td rowspan="3" style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:bold;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Raster</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Very good for representing remote sensing data</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Less aesthetic pleasing</span></p></td></tr><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Performs well with map algebra</span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Lacks attribute data flexibility<br></span></p></td></tr><tr><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;"></span></p></td><td style="width:288px;height:18px;background-color:transparent;vertical-align: middle;transform: rotate(0deg);border-bottom: 1px solid rgba(51, 51, 51, 1.00);border-top: 1px solid rgba(51, 51, 51, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);margin-bottom:0;margin-top:0;margin-left:0;margin-right:0;"><p style="margin:0;text-align:left;border-bottom: 0 solid rgba(0, 0, 0, 1.00);border-top: 0 solid rgba(0, 0, 0, 1.00);border-left: 0 solid rgba(0, 0, 0, 1.00);border-right: 0 solid rgba(0, 0, 0, 1.00);padding-bottom:2px;padding-top:2px;padding-left:5px;padding-right:5px;background-color:transparent;"><span style="font-family:'Arial';font-size:20px;font-weight:normal;font-style:normal;text-decoration:none;color:rgba(17, 17, 17, 1.00);background-color:transparent;">Data size increases with the increase in resolution</span></p></td></tr></tbody></table></div>
---
class: center
# Importing Spatial Data
.center[
.img-med[
]]
---
Now that we've learned about different types of spatial data, let's load some. We'll use the `rgdal` package for this:
```r
pacman::p_load(rgdal)
pkmng_countries <- readOGR(dsn = "Data",
layer = "pkmng_polygons")
```
```
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/marianovosolov/Dropbox/Studies/ecodatasci_tlv/spatial_analysis/Data", layer: "pkmng_polygons"
## with 246 features
## It has 12 fields
```
---
```r
head(pkmng_countries)
```
```
## class : SpatialPolygonsDataFrame
## features : 6
## extent : -61.88722, 50.37499, -18.01639, 42.61805 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 12
## # A tibble: 6 x 12
## FIPS ISO2 ISO3 UN NAME AREA POP2005 REGION SUBREGI LON
## <fct> <fct> <fct> <int> <fct> <int> <dbl> <int> <int> <dbl>
## 1 AC AG ATG 28 Anti~ 44 8.30e4 19 29 -61.8
## 2 AG DZ DZA 12 Alge~ 238174 3.29e7 2 15 2.63
## 3 AJ AZ AZE 31 Azer~ 8260 8.35e6 142 145 47.4
## 4 AL AL ALB 8 Alba~ 2740 3.15e6 150 39 20.1
## 5 AM AM ARM 51 Arme~ 2820 3.02e6 142 145 44.6
## 6 AO AO AGO 24 Ango~ 124670 1.61e7 2 17 17.5
## # ... with 2 more variables: LAT <dbl>, pkmn_rc <int>
```
---
As usual, there is a `tidyverse` alternative to view spatial objects and edit the associated data
.small[
```r
pacman::p_load(spdplyr)
pkmng_countries
```
```
## class : SpatialPolygonsDataFrame
## features : 246
## extent : -180, 180, -90, 83.57027 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 12
## # A tibble: 246 x 12
## FIPS ISO2 ISO3 UN NAME AREA POP2005 REGION SUBREGI LON
## <fct> <fct> <fct> <int> <fct> <int> <dbl> <int> <int> <dbl>
## 1 AC AG ATG 28 Anti… 44 8.30e4 19 29 -61.8
## 2 AG DZ DZA 12 Alge… 238174 3.29e7 2 15 2.63
## 3 AJ AZ AZE 31 Azer… 8260 8.35e6 142 145 47.4
## 4 AL AL ALB 8 Alba… 2740 3.15e6 150 39 20.1
## 5 AM AM ARM 51 Arme… 2820 3.02e6 142 145 44.6
## 6 AO AO AGO 24 Ango… 124670 1.61e7 2 17 17.5
## 7 AQ AS ASM 16 Amer… 20 6.41e4 9 61 -171.
## 8 AR AR ARG 32 Arge… 273669 3.87e7 19 5 -65.2
## 9 AS AU AUS 36 Aust… 768230 2.03e7 9 53 136.
## 10 BA BH BHR 48 Bahr… 71 7.25e5 142 145 50.6
## # … with 236 more rows, and 2 more variables: LAT <dbl>,
## # pkmn_rc <int>
```
]
---
With this package, you can use `tidyverse` functions to manipulate your spatial data:
.small[
```r
#filter out all countries without any Pokemon sightings, and let's rename the column while we're at it:
pkmng_countries %>%
filter(pkmn_rc > 0) %>%
rename(Country = NAME)
```
```
## class : SpatialPolygonsDataFrame
## features : 72
## extent : -179.999, 179.999, -55.9175, 83.10942 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 12
## # A tibble: 72 x 12
## FIPS ISO2 ISO3 UN Country AREA POP2005 REGION SUBREGI
## <fct> <fct> <fct> <int> <fct> <int> <dbl> <int> <int>
## 1 AO AO AGO 24 Angola 124670 1.61e7 2 17
## 2 AR AR ARG 32 Argent~ 273669 3.87e7 19 5
## 3 AS AU AUS 36 Austra~ 768230 2.03e7 9 53
## 4 BK BA BIH 70 Bosnia~ 5120 3.92e6 150 39
## 5 BM MM MMR 104 Burma 65755 4.80e7 142 35
## 6 BR BR BRA 76 Brazil 845942 1.87e8 19 5
## 7 BX BN BRN 96 Brunei~ 527 3.74e5 142 35
## 8 CA CA CAN 124 Canada 909351 3.23e7 19 21
## 9 CB KH KHM 116 Cambod~ 17652 1.40e7 142 35
## 10 CI CL CHL 152 Chile 74880 1.63e7 19 5
## # ... with 62 more rows, and 3 more variables: LON <dbl>, LAT <dbl>,
## # pkmn_rc <int>
```
]
---
Now let's load our point data:
.small[
```r
pkmng_points_dat <- read_csv("Data/pkmng_points.csv")
pkmng_points_dat
```
```
## # A tibble: 29,602 x 207
## pokemonId latitude longitude appearedLocalTime cellId_90m
## <dbl> <dbl> <dbl> <dttm> <dbl>
## 1 16 53.7 -0.443 2016-09-04 09:31:44 5.22e18
## 2 60 46.0 8.96 2016-09-06 07:57:06 5.15e18
## 3 19 46.2 15.3 2016-09-04 08:33:55 5.14e18
## 4 16 50.7 14.5 2016-09-04 09:23:24 5.12e18
## 5 23 44.5 11.3 2016-09-04 16:02:16 5.15e18
## 6 120 50.2 14.7 2016-09-03 15:51:17 5.12e18
## 7 46 33.6 -84.3 2016-09-07 01:05:11 9.87e18
## 8 10 46.0 8.95 2016-09-06 06:22:14 5.15e18
## 9 27 33.9 -118. 2016-09-06 05:13:44 9.28e18
## 10 23 34.3 -119. 2016-09-06 14:49:44 9.29e18
## # ... with 29,592 more rows, and 202 more variables:
## # cellId_180m <dbl>, cellId_370m <dbl>, cellId_730m <dbl>,
## # cellId_1460m <dbl>, cellId_2920m <dbl>, cellId_5850m <dbl>,
## # appearedTimeOfDay <chr>, appearedHour <dbl>,
## # appearedMinute <dbl>, appearedDayOfWeek <chr>, appearedDay <dbl>,
## # appearedMonth <dbl>, appearedYear <dbl>, terrainType <dbl>,
## # closeToWater <lgl>, city <chr>, continent <chr>, weather <chr>,
## # temperature <dbl>, windSpeed <dbl>, windBearing <dbl>,
## # pressure <dbl>, weatherIcon <chr>, sunriseMinutesMidnight <dbl>,
## # sunriseHour <dbl>, sunriseMinute <dbl>,
## # sunriseMinutesSince <dbl>, sunsetMinutesMidnight <dbl>,
## # sunsetHour <dbl>, sunsetMinute <dbl>, sunsetMinutesBefore <dbl>,
## # population_density <dbl>, urban <lgl>, suburban <lgl>,
## # midurban <lgl>, rural <lgl>, gymDistanceKm <dbl>,
## # gymIn100m <lgl>, gymIn250m <lgl>, gymIn500m <lgl>,
## # gymIn1000m <lgl>, gymIn2500m <lgl>, gymIn5000m <lgl>,
## # pokestopDistanceKm <dbl>, pokestopIn100m <lgl>,
## # pokestopIn250m <lgl>, pokestopIn500m <lgl>,
## # pokestopIn1000m <lgl>, pokestopIn2500m <lgl>,
## # pokestopIn5000m <lgl>, cooc_1 <lgl>, cooc_2 <lgl>, cooc_3 <lgl>,
## # cooc_4 <lgl>, cooc_5 <lgl>, cooc_6 <lgl>, cooc_7 <lgl>,
## # cooc_8 <lgl>, cooc_9 <lgl>, cooc_10 <lgl>, cooc_11 <lgl>,
## # cooc_12 <lgl>, cooc_13 <lgl>, cooc_14 <lgl>, cooc_15 <lgl>,
## # cooc_16 <lgl>, cooc_17 <lgl>, cooc_18 <lgl>, cooc_19 <lgl>,
## # cooc_20 <lgl>, cooc_21 <lgl>, cooc_22 <lgl>, cooc_23 <lgl>,
## # cooc_24 <lgl>, cooc_25 <lgl>, cooc_26 <lgl>, cooc_27 <lgl>,
## # cooc_28 <lgl>, cooc_29 <lgl>, cooc_30 <lgl>, cooc_31 <lgl>,
## # cooc_32 <lgl>, cooc_33 <lgl>, cooc_34 <lgl>, cooc_35 <lgl>,
## # cooc_36 <lgl>, cooc_37 <lgl>, cooc_38 <lgl>, cooc_39 <lgl>,
## # cooc_40 <lgl>, cooc_41 <lgl>, cooc_42 <lgl>, cooc_43 <lgl>,
## # cooc_44 <lgl>, cooc_45 <lgl>, cooc_46 <lgl>, cooc_47 <lgl>,
## # cooc_48 <lgl>, cooc_49 <lgl>, cooc_50 <lgl>, ...
```
]
---
We are going to create a shapefile from this table.
```r
pacman::p_load(sp)
xy <- tibble(longitude = pkmng_points_dat$longitude,
latitude = pkmng_points_dat$latitude)
pkmng_points <- SpatialPointsDataFrame(coords = xy,
data = pkmng_points_dat)
pkmng_points
```
```
## class : SpatialPointsDataFrame
## features : 29602
## extent : -158.0198, 175.6162, -42.9845, 68.43698 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## variables : 207
## names : pokemonId, latitude, longitude, appearedLocalTime, cellId_90m, cellId_180m, cellId_370m, cellId_730m, cellId_1460m, cellId_2920m, cellId_5850m, appearedTimeOfDay, appearedHour, appearedMinute, appearedDayOfWeek, ...
## min values : 1, 1.282526, 0.000046, 2016-09-02 21:50:11, 1.008879e+19, 1.008879e+19, 1.008879e+19, 1.008879e+19, 1.008879e+19, 1.008879e+19, 1.008879e+19, afternoon, 0, 0, dummy_day, ...
## max values : 149, -42.984504, -158.019758, 2016-09-08 03:51:27, 9.937698e+18, 9.937698e+18, 9.937698e+18, 9.937698e+18, 9.937698e+18, 9.937698e+18, 9.937698e+18, night, 23, 59, Wednesday, ...
```
---
# GeoJSON
So far we've been working with ESRI shapefiles. This is a format to encode geographic data, used widely in GIS systems.
Alternatively, there is the java-based GeoJSON format for encoding geographic data, which can support the features we've discussed (points, polygons).
There are advantages to GeoJSON, including reduced file sizes. Compare:
---
```r
pacman::p_load(geojsonio)
pacman::p_load(rmapshaper)
#we now convert our ESRI shapefile (a SpatialPoylgonsDataFrame object) to GeoJSON format:
pkmng_countries_json <- geojson_json(pkmng_countries) #this can be VERY time-consuming
```
Compare the file sizes: 6.4 Mb for the ESRI format, compared to 1.1 Mb for the GeoJSON format!
Read more about GeoJSON here: https://macwright.org/2015/03/23/geojson-second-bite.html
---
# Rasters
Finally, importing rasters is easily done with the `raster` package:
```r
pacman::p_load(raster)
pkmng_raster <- raster("Data/pkmng_raster.tif")
pkmng_raster
```
```
## class : RasterLayer
## dimensions : 87, 180, 15660 (nrow, ncol, ncell)
## resolution : 2, 2 (x, y)
## extent : -180, 180, -90.42973, 83.57027 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /Users/marianovosolov/Dropbox/Studies/ecodatasci_tlv/spatial_analysis/Data/pkmng_raster.tif
## names : pkmng_raster
## values : 0, 9.535896 (min, max)
```
---
# Projections
One of the most important things to remember when dealing with spatial data is the projection.
A projection is a mathematical transformation of coordinates from the surface of a sphere on a 2D plane.
.center[]
---
# Types of projections
## Projections by surface
Classification based on the type of surface onto which the globe is projected:
* Cylindrical
* Conic
* Azimuthal
etc.
---
## Projections by preservation of a metric property
Classification based on which property is maintained in the projection:
* Conformal (preserves shape, usually distorts area)
* Equal-area (preserves area, usually distorts shape)
* Equidistant (preserves distance from some point or line)
etc.
---
You can view the Coordinate Reference System of your spatial objects:
```r
crs(pkmng_countries)
```
```
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
```
You can see different values here, including the datum, units, etc.
---
Notice our points shapefile has no projection, because we didn't set it when creating the object:
```r
crs(pkmng_points)
```
```
## CRS arguments: NA
```
---
We can set the projection for this object, either by writing a string of our desired projection, or extracting the projection from an existing object:
```r
proj1 <- crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj1
```
```
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
```
```r
proj2 <- crs(pkmng_countries)
proj2
```
```
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
```
---
```r
proj4string(pkmng_points) <- proj1
pkmng_points
```
```
## class : SpatialPointsDataFrame
## features : 29602
## extent : -158.0198, 175.6162, -42.9845, 68.43698 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 207
## names : pokemonId, latitude, longitude, appearedLocalTime, cellId_90m, cellId_180m, cellId_370m, cellId_730m, cellId_1460m, cellId_2920m, cellId_5850m, appearedTimeOfDay, appearedHour, appearedMinute, appearedDayOfWeek, ...
## min values : 1, -42.984504, -158.019758, 1472853011, 43172735374852096, 43172735576178688, 43172735307743232, 43172732086517760, 43172719201615872, 43172667662008320, 43172598942531584, afternoon, 0, 0, dummy_day, ...
## max values : 149, 68.43698, 175.616223, 1473306687, 12280882455157669888, 12280882455224778752, 12280882455493214208, 1.2280882456567e+19, 12280882443682054144, 12280882495221661696, 12280882701380091904, night, 23, 59, Wednesday, ...
```
---
You can also transform already existing objects to have different projections:
```r
pkmng_points_p <- spTransform(pkmng_points,
crs("+proj=laea +lat_0=0 +lon_0=0"))
pkmng_points_p
```
```
## class : SpatialPointsDataFrame
## features : 29602
## extent : -9285067, 11178209, -11993372, 9520128 (xmin, xmax, ymin, ymax)
## crs : +proj=laea +lat_0=0 +lon_0=0 +ellps=WGS84
## variables : 207
## names : pokemonId, latitude, longitude, appearedLocalTime, cellId_90m, cellId_180m, cellId_370m, cellId_730m, cellId_1460m, cellId_2920m, cellId_5850m, appearedTimeOfDay, appearedHour, appearedMinute, appearedDayOfWeek, ...
## min values : 1, -42.984504, -158.019758, 1472853011, 43172735374852096, 43172735576178688, 43172735307743232, 43172732086517760, 43172719201615872, 43172667662008320, 43172598942531584, afternoon, 0, 0, dummy_day, ...
## max values : 149, 68.43698, 175.616223, 1473306687, 12280882455157669888, 12280882455224778752, 12280882455493214208, 1.2280882456567e+19, 12280882443682054144, 12280882495221661696, 12280882701380091904, night, 23, 59, Wednesday, ...
```
---
Why are projections so important?
Here's our points mapped on top of a shapefile of countries:
```r
plot(countries, axes = T)
plot(pkmng_points, add = T, pch = 1, col = "red")
```
<img src="index_files/figure-html/map-1.png" style="display: block; margin: auto;" />
---
Now, when the projections do not match:
```r
countries_p <- spTransform(countries,
crs("+proj=laea +lat_0=0 +lon_0=0"))
plot(countries_p, axes = T)
plot(pkmng_points, add = T, pch = 1, col = "red")
```
<img src="index_files/figure-html/map 2-1.png" style="display: block; margin: auto;" />
---
Whereas matching projections give us all the coordinates in the right place:
```r
plot(countries_p, axes = T)
plot(pkmng_points_p, add = T, pch = 1, col = "red")
```
<img src="index_files/figure-html/map 3-1.png" style="display: block; margin: auto;" />
---
# Extent
Spatial objects also have extents - these are the minimum and maximum X and Y coordinates of the object. They don't have to have data in them. The coordinates of course depend on your projection and units:
```r
countries@bbox
```
```
## min max
## x -180 180.00000
## y -90 83.57027
```
```r
countries_p@bbox
```
```
## min max
## x -12726751 12733710
## y -12641498 12299760
```
---
For raster objects, you find the extent with this syntax:
```r
pkmng_raster@extent
```
```
## class : Extent
## xmin : -180
## xmax : 180
## ymin : -86.42973
## ymax : 83.57027
```
---
You can change the extent of objects using the `crop()` or `extend()` functions, depending on whether you want to decrease or increase the extent:
```r
pkmng_raster_c <- crop(pkmng_raster,
extent(-50,
100,
0,
80))
```
---
```r
plot(pkmng_raster)
plot(pkmng_raster_c, add = T, col = "red")
```
<img src="index_files/figure-html/extent 4-1.png" style="display: block; margin: auto;" />
All the red cells are in the new, limited extent!
---
## Why is this important?
If you try to do spatial arithmetics on objects with mismatching extents and projections, you'll get errors. For example:
```r
raster::intersect(countries_p, pkmng_raster_c)
```
```
## Warning in intersect(x, y): non identical CRS
```
```
## Warning in intersect(x, y): polygons do not intersect
```
```
## NULL
```
---
class: center, middle
# Data manipulation
---
class: center, middle
# Geoprocessing
a GIS operation used to manipulate spatial data. A typical geoprocessing operation takes an input dataset, performs an operation on that dataset, and returns the result of the operation as an output dataset
---
### **Buffer points** - creating a polygon based on proximity
**For example:** buffer some of the pokemon points to 100km radius
We'll first subset the data to have less points
```r
sub_pkmng_p<- subset(pkmng_points,pokemonId %in%c(8:10))
plot(countries)
plot(sub_pkmng_p,add = T, col = "red")
```
<img src="index_files/figure-html/unnamed-chunk-2-1.png" style="display: block; margin: auto;" />
---
Now we'll buffer the points in a 100km radius to make them into polygons
```r
buffer_pkmng_points<- raster::buffer(sub_pkmng_p,width = 100000)
plot(countries)
plot(buffer_pkmng_points,add = T, col = "red")
```
<img src="index_files/figure-html/unnamed-chunk-3-1.png" style="display: block; margin: auto;" />
---
### **Clip** - an overlay function that cuts out an input layer with the extent of a defined feature boundary. The result of this tool is a new clipped output layer.
There are a few ways of clipping polygons in R, we will look at using the rgdal package.
---
Lets clip the buffered points to to have the ones that fall in Italy and to fit the boundaries of the Italian coastline.
First, well filter the countries shapefile to have a shapefile of Italy
```r
Italy <- countries %>%
filter(NAME == "Italy")
plot(Italy)
plot(buffer_pkmng_points,add = T, col = "red")
```
<img src="index_files/figure-html/unnamed-chunk-4-1.png" style="display: block; margin: auto;" />
---
# `rgeos` package example
Now, lets clip the the buffered points based on the Italian shapefile we created
```r
clip_pkmng_rgeos<- rgeos::gIntersection(Italy, buffer_pkmng_points, byid = TRUE, drop_lower_td = TRUE)
plot(Italy)
plot(clip_pkmng_rgeos,add = T, col = "red")
```
<img src="index_files/figure-html/unnamed-chunk-5-1.png" style="display: block; margin: auto;" />
---
### **Intersect** - similar to the clip tool because the extents of input features defines the output. The only exception is that it preserves attributes from all the data sets that overlap each other in the output
Lets intersect the points with the country shapefile to find which points belong to which country
```r
intersect_pkmg<- raster::intersect(pkmng_points,countries)
```
---
### **Merge** - combines data sets that are the same data type
```r
#subset the countries
sub_count<- countries %>% filter(NAME==c("France","Italy"))
sub_count_sf<- sf::st_as_sf(sub_count)
merge_count<- sf::st_combine(sub_count_sf)
plot(merge_count)
```
<img src="index_files/figure-html/unnamed-chunk-7-1.png" style="display: block; margin: auto;" />
---
### **Dissolve** - The Dissolve Tool unifies boundaries based on common attribute values
```r
diss_sub_count<- rgeos::gUnaryUnion(sub_count)
plot(diss_sub_count)
```
<img src="index_files/figure-html/unnamed-chunk-8-1.png" style="display: block; margin: auto;" />
---
# Zonal statistics
Vector-Raster analysis. Allows you to calculate statistics from raster data for each feature in the vector data
**For Example:** We want to calculate the temperature in each point that a pokemon was collected
---
Download and plot the mean annual temperature from the web using the getData function
```r
climate <- raster::getData('worldclim', var='bio', res=10)
annual_prec<- climate$bio12
plot(annual_prec)
```
<img src="index_files/figure-html/unnamed-chunk-9-1.png" style="display: block; margin: auto;" />
---
Now lets calculate the mean temperature in each point
```r
poke_annual_prec<- raster::extract(annual_prec, pkmng_points,fun = mean,na.rm=T,df=T)
head(poke_annual_prec)
```
```
## ID bio12
## 1 1 663
## 2 2 1326
## 3 3 1095
## 4 4 585
## 5 5 792
## 6 6 492
```
---
The data has no identification of which point it belongs to. To have that we need to merge it with the points identifiers
```r
poke_annual_prec_final<- cbind(pkmng_points@data$pokemonId,poke_annual_prec$bio12)
names(poke_annual_prec_final)<- c("pokemonId","annual_prec")
head(poke_annual_prec_final)
```
```
## [,1] [,2]
## [1,] 16 663
## [2,] 60 1326
## [3,] 19 1095
## [4,] 16 585
## [5,] 23 792
## [6,] 120 492
```
---
# Spatial autocorrelation
Spatial data are not independent - violating assumptions of statistical tests.
Spatial autocorrelation decreases with distance - but the distance depends on your data.
So how do we test for spatial autocorrelation and what can we do about it?
---
## Testing for spatial autocorrelation
We'll use the `pgirmess` package:
```r
pacman::p_load(pgirmess)
#first we'll take a small sample of the points data, because this test is very time-consuming
pkmng_points_dat_s <- sample_frac(pkmng_points_dat, 0.05)
pgi.cor <- correlog(coords = pkmng_points_dat_s[,c("longitude","latitude")], #longitude and latitude, in that order
z = pkmng_points_dat_s$temperature, #the variable we want to test for autocorrelation in
method = "Moran")
```
---
The X axis is distance, and the Y axis is our measure of autocorrelation. Values higher than 0 mean there is spatial autocorrelation in that distance class.
```r
plot(pgi.cor)
```
<img src="index_files/figure-html/sp ac 2-1.png" style="display: block; margin: auto;" />
---
And this is what it looks like when there's no spatial autocorrelation:
```r
pgi.cor2 <- correlog(coords = pkmng_points_dat_s[,c("longitude","latitude")], #longitude and latitude, in that order
z = pkmng_points_dat_s$pokemonId, #the variable we want to test for autocorrelation in
method = "Moran")
plot(pgi.cor2)
```
<img src="index_files/figure-html/sp ac 3-1.png" style="display: block; margin: auto;" />
---
## Dealing with spatial autocorrelation
There are various methods to "eliminate" spatial autocorrelation for hypothesis testing and statistical inference. All are flawed, each in their own way.
We won't talk about them here, but you can read more in the literature:
https://onlinelibrary.wiley.com/doi/10.1111/j.2007.0906-7590.05171.x
https://onlinelibrary.wiley.com/doi/10.1111/jbi.12953
---
# Spatial non-stationarity
Often, relationships between variables can shift in space. So global coefficients of models can be very poor at describing the actual behavior of your data.
For example, let's look for a relationship between wind speed and temperature in our Pokemon Go sightings data:
.small[
```r
lm.model <- lm(windSpeed ~
temperature,
data = pkmng_points_dat_s)
summary(lm.model)
```
```
##
## Call:
## lm(formula = windSpeed ~ temperature, data = pkmng_points_dat_s)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.071 -3.112 -0.737 2.241 49.873
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.06517 0.40037 12.651 < 2e-16 ***
## temperature 0.07405 0.01766 4.194 2.91e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.249 on 1478 degrees of freedom
## Multiple R-squared: 0.01176, Adjusted R-squared: 0.01109
## F-statistic: 17.59 on 1 and 1478 DF, p-value: 2.906e-05
```
]
---
Now, we'll map the residuals from this model:
```r
resids<-residuals(lm.model)
colours <- c("dark blue", "blue", "red", "dark red")
map.resids <- SpatialPointsDataFrame(data = data.frame(resids),
coords = cbind(pkmng_points_dat_s$longitude,
pkmng_points_dat_s$latitude)
)
```
---
```r
spplot(map.resids, cuts=quantile(resids), col.regions=colours, cex=1)
```
<img src="index_files/figure-html/gwr 3-1.png" style="display: block; margin: auto;" />
Notice how the red and blue dots are clustered together, and not randomly distributed - showing a spatial pattern in the residuals.
---
We can use a method called GWR (Globally Weighted Regression) to account for this - basically, you calculate a regression coefficient for each datum, weighing the contribution of all data based on their spatial distance from your focal datum.
.small[
```r
pacman::p_load(spgwr)
GWRbandwidth <- gwr.sel(windSpeed ~
temperature,
data = pkmng_points_dat_s,
coords = cbind(pkmng_points_dat_s$longitude,
pkmng_points_dat_s$latitude),
adapt = T)
```
```
## Adaptive q: 0.381966 CV score: 25971.6
## Adaptive q: 0.618034 CV score: 26470.16
## Adaptive q: 0.236068 CV score: 25779.09
## Adaptive q: 0.145898 CV score: 25429.72
## Adaptive q: 0.09016994 CV score: 24672.23
## Adaptive q: 0.05572809 CV score: 23923.28
## Adaptive q: 0.03444185 CV score: 22690.01
## Adaptive q: 0.02128624 CV score: 21834.46
## Adaptive q: 0.01315562 CV score: 20909.26
## Adaptive q: 0.008130619 CV score: 20374.53
## Adaptive q: 0.005024999 CV score: 21990.94
## Adaptive q: 0.01005 CV score: 20553.34
## Adaptive q: 0.006944377 CV score: 21538.02
## Adaptive q: 0.008863756 CV score: 20508.88
## Adaptive q: 0.007677515 CV score: 20717.15
## Adaptive q: 0.008410652 CV score: 20342.56
## Adaptive q: 0.008357601 CV score: 20328.39
## Adaptive q: 0.008304636 CV score: 20328.57
## Adaptive q: 0.008357601 CV score: 20328.39
```
]
---
```r
gwr.model <- gwr(windSpeed ~
temperature,
data = pkmng_points_dat_s,
coords = cbind(pkmng_points_dat_s$longitude,
pkmng_points_dat_s$latitude),
adapt = GWRbandwidth,
hatmatrix = T,
se.fit = T
)
gwr.results <- as_tibble(gwr.model$SDF)
```
---
Now we can view the table of the results, and see the coefficient for each datum in our dataset - i.e. for each point locality. We can see they vary quite a bit.
.small[
```r
gwr.results
```
```
## # A tibble: 1,480 x 14
## sum.w X.Intercept. temperature X.Intercept._se temperature_se
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 30.7 31.2 -1.02 3.18 0.131
## 2 13.8 6.22 -0.105 6.27 0.257
## 3 14.4 11.4 -0.0643 7.80 0.390
## 4 13.2 2.73 0.317 3.33 0.246
## 5 63.8 5.35 0.0547 1.15 0.0454
## 6 45.0 6.98 0.0553 1.41 0.0763
## 7 20.7 10.4 -0.201 2.64 0.135
## 8 13.9 1.26 0.365 3.11 0.182
## 9 14.4 11.4 -0.0654 7.78 0.389
## 10 13.1 -9.07 0.666 3.64 0.158
## # ... with 1,470 more rows, and 9 more variables: gwr.e <dbl>,
## # pred <dbl>, pred.se <dbl>, localR2 <dbl>,
## # X.Intercept._se_EDF <dbl>, temperature_se_EDF <dbl>,
## # pred.se.1 <dbl>, coord.x <dbl>, coord.y <dbl>
```
]
---
Now let's plot this and see exactly what this non-stationarity looks like:
.small[
```r
pkmng_points_dat_s$coef <- gwr.results$temperature
countries_outline <- fortify(countries,
region = "NAME")
mapWorld <- borders("world", colour=NA, fill="gray")
```
```
## Warning: package 'maps' was built under R version 3.5.2
```
```r
gwr.map <- ggplot(data = pkmng_points_dat_s,
aes(x = longitude,
y = latitude,
colour = coef)
) +
mapWorld +
geom_point() +
scale_colour_gradient2(low = "red",
mid = "white",
high = "blue",
midpoint = 0,
space = "rgb",
na.value = "grey50",
guide = "colourbar",
guide_legend(title="Coefs")) +
coord_equal()
```
```
## Warning: Non Lab interpolation is deprecated
```
]
---
And see how the regression coefficients change with space - blue for positive coefficients and red for negative ones:
```r
gwr.map
```