forked from AdamWilsonLabEDU/SpatialDataScience
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy path04_Spatial_with_sf.html
More file actions
1073 lines (994 loc) · 50.3 KB
/
04_Spatial_with_sf.html
File metadata and controls
1073 lines (994 loc) · 50.3 KB
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">
<head>
<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />
<title>Working with Spatial Data</title>
<script src="site_libs/jquery-1.11.3/jquery.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/cosmo.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<script src="site_libs/jqueryui-1.11.4/jquery-ui.min.js"></script>
<link href="site_libs/tocify-1.9.1/jquery.tocify.css" rel="stylesheet" />
<script src="site_libs/tocify-1.9.1/jquery.tocify.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<link href="site_libs/highlightjs-9.12.0/textmate.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/highlight.js"></script>
<link href="site_libs/font-awesome-5.0.13/css/fa-svg-with-js.css" rel="stylesheet" />
<script src="site_libs/font-awesome-5.0.13/js/fontawesome-all.min.js"></script>
<script src="site_libs/font-awesome-5.0.13/js/fa-v4-shims.min.js"></script>
<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
pre:not([class]) {
background-color: white;
}
</style>
<script type="text/javascript">
if (window.hljs) {
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
if (document.readyState && document.readyState === "complete") {
window.setTimeout(function() { hljs.initHighlighting(); }, 0);
}
}
</script>
<style type="text/css">
h1 {
font-size: 34px;
}
h1.title {
font-size: 38px;
}
h2 {
font-size: 30px;
}
h3 {
font-size: 24px;
}
h4 {
font-size: 18px;
}
h5 {
font-size: 16px;
}
h6 {
font-size: 12px;
}
.table th:not([align]) {
text-align: left;
}
</style>
<link rel="stylesheet" href="styles.css" type="text/css" />
</head>
<body>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
code {
color: inherit;
background-color: rgba(0, 0, 0, 0.04);
}
img {
max-width:100%;
height: auto;
}
.tabbed-pane {
padding-top: 12px;
}
.html-widget {
margin-bottom: 20px;
}
button.code-folding-btn:focus {
outline: none;
}
</style>
<style type="text/css">
/* padding for bootstrap navbar */
body {
padding-top: 51px;
padding-bottom: 40px;
}
/* offset scroll position for anchor links (for fixed navbar) */
.section h1 {
padding-top: 56px;
margin-top: -56px;
}
.section h2 {
padding-top: 56px;
margin-top: -56px;
}
.section h3 {
padding-top: 56px;
margin-top: -56px;
}
.section h4 {
padding-top: 56px;
margin-top: -56px;
}
.section h5 {
padding-top: 56px;
margin-top: -56px;
}
.section h6 {
padding-top: 56px;
margin-top: -56px;
}
</style>
<script>
// manage active state of menu based on current page
$(document).ready(function () {
// active menu anchor
href = window.location.pathname
href = href.substr(href.lastIndexOf('/') + 1)
if (href === "")
href = "index.html";
var menuAnchor = $('a[href="' + href + '"]');
// mark it active
menuAnchor.parent().addClass('active');
// if it's got a parent navbar menu mark it active as well
menuAnchor.closest('li.dropdown').addClass('active');
});
</script>
<div class="container-fluid main-container">
<!-- tabsets -->
<script>
$(document).ready(function () {
window.buildTabsets("TOC");
});
</script>
<!-- code folding -->
<script>
$(document).ready(function () {
// move toc-ignore selectors from section div to header
$('div.section.toc-ignore')
.removeClass('toc-ignore')
.children('h1,h2,h3,h4,h5').addClass('toc-ignore');
// establish options
var options = {
selectors: "h1,h2",
theme: "bootstrap3",
context: '.toc-content',
hashGenerator: function (text) {
return text.replace(/[.\\/?&!#<>]/g, '').replace(/\s/g, '_').toLowerCase();
},
ignoreSelector: ".toc-ignore",
scrollTo: 0
};
options.showAndHide = true;
options.smoothScroll = true;
// tocify
var toc = $("#TOC").tocify(options).data("toc-tocify");
});
</script>
<style type="text/css">
#TOC {
margin: 25px 0px 20px 0px;
}
@media (max-width: 768px) {
#TOC {
position: relative;
width: 100%;
}
}
.toc-content {
padding-left: 30px;
padding-right: 40px;
}
div.main-container {
max-width: 1200px;
}
div.tocify {
width: 20%;
max-width: 260px;
max-height: 85%;
}
@media (min-width: 768px) and (max-width: 991px) {
div.tocify {
width: 25%;
}
}
@media (max-width: 767px) {
div.tocify {
width: 100%;
max-width: none;
}
}
.tocify ul, .tocify li {
line-height: 20px;
}
.tocify-subheader .tocify-item {
font-size: 0.90em;
padding-left: 25px;
text-indent: 0;
}
.tocify .list-group-item {
border-radius: 0px;
}
</style>
<!-- setup 3col/9col grid for toc_float and main content -->
<div class="row-fluid">
<div class="col-xs-12 col-sm-4 col-md-3">
<div id="TOC" class="tocify">
</div>
</div>
<div class="toc-content col-xs-12 col-sm-8 col-md-9">
<div class="navbar navbar-inverse navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="index.html">GEO 503: Spatial Data Science</a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
<li>
<a href="index.html">Home</a>
</li>
<li>
<a href="Syllabus.html">Syllabus</a>
</li>
<li>
<a href="Schedule.html">Schedule</a>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
Content
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="CourseContent.html">About Course Content</a>
</li>
<li class="dropdown-header">Module 1: Introduction to R</li>
<li>
<a href="00_CourseIntroductionFrame.html">00 Course Introduction</a>
</li>
<li>
<a href="01_Rintro.html">01 First Steps</a>
</li>
<li>
<a href="02_graphics.html">02 Graphics</a>
</li>
<li>
<a href="03_DataWrangling.html">03 Data Wrangling</a>
</li>
<li>
<a href="03b_DataWrangling.html">03 Data Wrangling 2</a>
</li>
<li class="divider"></li>
<li class="dropdown-header">Module 2: Spatial Analyses</li>
<li>
<a href="04_Spatial_with_sf.html">04 Spatial Data with sf</a>
</li>
<li>
<a href="05_Raster.html">05 Spatial Raster Data</a>
</li>
<li class="divider"></li>
<li class="dropdown-header">Module 3: Advanced Programming</li>
<li>
<a href="06_CreatingWorkflows.html">06 Creating Workflows</a>
</li>
<li>
<a href="07_Reproducibile.html">07 Reproducible Research</a>
</li>
<li>
<a href="08_WeatherData.html">08 Weather Data Processing</a>
</li>
<li>
<a href="09_RemoteSensing_appeears.html">09 Satellite Data Processing</a>
</li>
<li>
<a href="11_ParallelProcessing.html">11 Parallel Processing</a>
</li>
<li>
<a href="12_DynamicVisualization.html">12 Dynamic Visualization</a>
</li>
<li>
<a href="13_SDM_Exercise.html">13 Species Distribution Modeling</a>
</li>
<li class="divider"></li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
Assignments
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="Tasklist.html">Task list</a>
</li>
<li>
<a href="DataCamp.html">DataCamp</a>
</li>
<li>
<a href="PackageIntro.html">Package Introduction</a>
</li>
<li>
<a href="Project.html">Final Project</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
Resources
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="Provenance.html">Provenance</a>
</li>
<li>
<a href="Projects_2018.html">2018 Project Drafts</a>
</li>
<li>
<a href="Projects_2017.html">2017 Final Projects</a>
</li>
<li>
<a href="GitSSHNotes.html">Setting up Github</a>
</li>
</ul>
</li>
</ul>
<ul class="nav navbar-nav navbar-right">
<li>
<a href="https://github.com/adammwilson/RDataScience">
<span class="fa fa-github fa-lg"></span>
</a>
</li>
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<div class="fluid-row" id="header">
<h1 class="title toc-ignore">Working with Spatial Data</h1>
</div>
<div class="extraswell">
<button data-toggle="collapse" class="btn btn-link" data-target="#pres">
View Presentation
</button>
<a href="presentations/PS_08_sf.html" target="_blank">Open presentation in a new tab</a>
<div id="pres" class="collapse">
<div class="embed-responsive embed-responsive-16by9">
<p><iframe class="embed-responsive-item" src="presentations/PS_08_sf.html" allowfullscreen></iframe> <em>Click on presentation and then use the space bar to advance to the next slide or escape key to show an overview.</em></p>
</div>
</div>
</div>
<div id="download" class="section level2">
<h2>Download</h2>
<table>
<thead>
<tr class="header">
<th align="center"><a href="scripts/04_Spatial_with_sf_nocomments.R"><i class="fas fa-code fa-2x" aria-hidden="true"></i><br> R Script</a></th>
<th align="center"><a href="scripts/04_Spatial_with_sf.R"><i class="fa fa-file-code-o fa-2x"></i> <br> Commented R Script</a></th>
<th align="center"><a href="scripts/04_Spatial_with_sf.Rmd"><i class="far fa-file-alt fa-2x"></i> <br> Rmd Script</a></th>
</tr>
</thead>
<tbody>
</tbody>
</table>
</div>
<div id="spatial-packages" class="section level1">
<h1>Spatial packages</h1>
<p>In R, there are two main lineages of tools for dealing with spatial data: sp and sf.</p>
<ul>
<li><p><code>sp</code> has been around for a while (the first release was in 2005), and it has a rich ecosystem of tools built on top of it. However, it uses a rather complex data structure, which can make it challenging to use.</p></li>
<li><p><code>sf</code> is newer (first released in October 2016!) so it doesn’t have such a rich ecosystem. However, it’s much easier to use and fits in very naturally with the tidyverse, and the ecosystem around it will grow rapidly.</p></li>
</ul>
<p>In this class, we’re going to use sf, so start by installing it(if you haven’t already):</p>
<pre class="r"><code>install.packages("sf")</code></pre>
<pre class="r"><code>library(tidyverse)
library(sf)
library(maps)</code></pre>
<p><a href="https://en.wikipedia.org/wiki/Simple_Features">Simple features</a> or <a href="http://www.opengeospatial.org/standards/sfa"><em>simple feature access</em></a> refers to a formal standard (ISO 19125-1:2004) that describes how objects in the real world can be represented in computers, with emphasis on the <em>spatial</em> geometry of these objects. It also describes how such objects can be stored in and retrieved from databases, and which geometrical operations should be defined for them.</p>
<p>The standard is widely implemented in spatial databases (such as PostGIS), commercial GIS (e.g., <a href="http://www.esri.com/">ESRI ArcGIS</a>) and forms the vector data basis for libraries such as <a href="http://www.gdal.org/">GDAL</a>. A subset of simple features forms the <a href="http://geojson.org/">GeoJSON</a> standard.</p>
<p>If you work with PostGis or GeoJSON you may have come across the <a href="https://en.wikipedia.org/wiki/Well-known_text">WKT (well-known text)</a> format, for example like these:</p>
<pre><code>POINT (30 10)
LINESTRING (30 10, 10 30, 40 40)
POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))</code></pre>
<p><code>sf</code> implements this standard natively in R. Data are structured and conceptualized very differently from the <code>sp</code> approach.</p>
<p>In <code>sf</code> spatial objects are stored as a simple data frame with a special column that contains the information for the geographic coordinates. That special column is a list with the same length as the number of rows in the data frame. Each of the individual list elements then can be of any length needed to hold the coordinates that correspond to an individual feature.</p>
</div>
<div id="data-io" class="section level1">
<h1>Data I/O</h1>
<div id="loading-data" class="section level2">
<h2>Loading data</h2>
<p>To read spatial data in R, use <code>read_sf()</code>. The following example reads an example dataset built into the sf package:</p>
<pre class="r"><code># The counties of North Carolina
file=system.file("shape/nc.shp", package = "sf")
file</code></pre>
<pre><code>## [1] "/Library/Frameworks/R.framework/Versions/3.5/Resources/library/sf/shape/nc.shp"</code></pre>
<pre class="r"><code>nc <- read_sf(file,
quiet = T,
stringsAsFactors = FALSE
)</code></pre>
<p>I recommend always setting <code>quiet = TRUE</code> and <code>stringsAsFactors = FALSE</code>.</p>
<p>Here we’re loading from a <strong>shapefile</strong> which is the way spatial data is most commonly stored. Despite the name a shapefile isn’t just one file, but is a collection of files that have the same name, but different extensions. Typically you’ll have four files:</p>
<ul>
<li><code>.shp</code> contains the geometry, and <code>.shx</code> contains an index into that geometry.</li>
<li><code>.dbf</code> contains metadata about each geometry (the other columns in the data frame).</li>
<li><code>.prf</code> contains the coordinate system and projection information. You’ll learn more about that shortly.</li>
</ul>
<p><code>read_sf()</code> can read in the majority of spatial file formats, so don’t worry if your data isn’t in a shapefile; the chances are <code>read_sf()</code> will still be able to read it.</p>
</div>
<div id="converting-data" class="section level2">
<h2>Converting data</h2>
<p>If you get a spatial object created by another package, us <code>st_as_sf()</code> to convert it to sf. For example, you can take data from the maps package (included in base R) and convert it to sf:</p>
<pre class="r"><code>nz_map <- maps::map("nz", plot = FALSE) #load a 'map' in non-sf format
nz_sf <- st_as_sf(nz_map) # convert it to sf format</code></pre>
</div>
</div>
<div id="data-structure" class="section level1">
<h1>Data structure</h1>
<p><code>nc</code> is a data frame, and not a tibble, so when printing, it’s a good idea to use <code>head()</code> so you only see the first few rows:</p>
<pre class="r"><code>head(nc)</code></pre>
<pre><code>## Simple feature collection with 6 features and 14 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
## epsg (SRID): 4267
## proj4string: +proj=longlat +datum=NAD27 +no_defs
## # A tibble: 6 x 15
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## <dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <dbl>
## 1 0.114 1.44 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.23 1827 1827 Alle… 37005 37005 3 487 0
## 3 0.143 1.63 1828 1828 Surry 37171 37171 86 3188 5
## 4 0.07 2.97 1831 1831 Curr… 37053 37053 27 508 1
## 5 0.153 2.21 1832 1832 Nort… 37131 37131 66 1421 9
## 6 0.097 1.67 1833 1833 Hert… 37091 37091 46 1452 7
## # ... with 5 more variables: NWBIR74 <dbl>, BIR79 <dbl>, SID79 <dbl>,
## # NWBIR79 <dbl>, geometry <MULTIPOLYGON [°]></code></pre>
<pre class="r"><code>head(nz_sf)</code></pre>
<pre><code>## Simple feature collection with 6 features and 1 field
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 166.458 ymin: -46.91705 xmax: 175.552 ymax: -36.09273
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## geometry ID
## 1 POLYGON ((166.458 -45.93695... Anchor.Island
## 2 POLYGON ((174.2599 -41.2092... Arapawa.Island
## 3 POLYGON ((166.58 -46.31315,... Coal.Island
## 4 POLYGON ((167.5798 -46.8738... Codfish.Island
## 5 POLYGON ((173.9064 -40.8492... D'Urville.Island
## 6 POLYGON ((175.5359 -36.3915... Great.Barrier.Island</code></pre>
<p>This is an ordinary data frame, with one exception: the <strong>geometry</strong> column. This column contains <strong>simple features</strong>, a standard way of representing two dimesional geometries like points, lines, polygons, multilines, and multipolygons. Multilines and multipolygons are nededed to represent geographic phenomena like a river with multiple branches, or a state made up of multiple islands.</p>
<pre class="r"><code>nc$geometry</code></pre>
<pre><code>## Geometry set for 100 features
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID): 4267
## proj4string: +proj=longlat +datum=NAD27 +no_defs
## First 5 geometries:</code></pre>
<pre><code>## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...</code></pre>
<pre><code>## MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...</code></pre>
<pre><code>## MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...</code></pre>
<pre><code>## MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...</code></pre>
<pre><code>## MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...</code></pre>
<p>Use <code>plot()</code> to show the geometry. You’ll learn how to use ggplot2 for more complex data visualisations later.</p>
<pre class="r"><code>plot(nc$geometry)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-8-1.png" /><!-- --></p>
<div id="manipulating-with-dplyr" class="section level2">
<h2>Manipulating with dplyr</h2>
<p>Since an <code>sf</code> object is just a data frame, you can manipulate it with dplyr. The following example gives you a taste:</p>
<pre class="r"><code>nz_sf %>%
mutate(area = as.numeric(st_area(geometry))) %>%
filter(area > 1e10)</code></pre>
<pre><code>## Simple feature collection with 2 features and 2 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 166.3961 ymin: -46.74155 xmax: 178.5629 ymax: -34.39895
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## ID area geometry
## 1 North.Island 113469632351 POLYGON ((172.7433 -34.4421...
## 2 South.Island 150444467051 POLYGON ((172.6391 -40.5135...</code></pre>
<p><code>st_area()</code> returns an object with units (i.e. <em>m</em><sup>2</sup>), which is precise, but a little annoying to work with. I used <code>as.numeric()</code> to convert to a regular numeric vector.</p>
</div>
<div id="geometry" class="section level2">
<h2>Geometry</h2>
<p>The geometry column is a list-column. In brief, they’re the richest and most complex type of column because a list can contain any other data structure, including other lists.</p>
<p>It’s worthwhile to pull out one piece so you can see what’s going on under the hood:</p>
<pre class="r"><code>str(nc$geometry[[1]])</code></pre>
<pre><code>## List of 1
## $ :List of 1
## ..$ : num [1:27, 1:2] -81.5 -81.5 -81.6 -81.6 -81.7 ...
## - attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"</code></pre>
<pre class="r"><code>plot(nc$geometry[[1]])</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-10-1.png" /><!-- --></p>
<p>Note the use of <code>[[</code> to extract a single element, here, the first polygon.</p>
<p>This is a list of lists of matrices:</p>
<ul>
<li>The top-level list has one element for each “landmass” in the county. However, we can find a more interesting case:</li>
</ul>
<pre class="r"><code>n <- nc$geometry %>% map_int(length)
table(n)</code></pre>
<pre><code>## n
## 1 2 3
## 94 4 2</code></pre>
<pre class="r"><code> interesting <- nc$geometry[n == 3][[1]]
plot(interesting)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-11-1.png" /><!-- --></p>
<pre class="r"><code>str(interesting)</code></pre>
<pre><code>## List of 3
## $ :List of 1
## ..$ : num [1:26, 1:2] -76 -76 -76 -76 -76.1 ...
## $ :List of 1
## ..$ : num [1:7, 1:2] -76 -76 -75.9 -75.9 -76 ...
## $ :List of 1
## ..$ : num [1:5, 1:2] -75.9 -75.9 -75.8 -75.8 -75.9 ...
## - attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"</code></pre>
<ul>
<li>This is a county made up of three non-contiguous pieces.</li>
<li>The second-level list is not used in this dataset, but is needed when you have a landmass that contains a lake. (Or a landmass that contains a lake which has an island which has a pond).</li>
<li>Each row of the matrix gives the location of a point on the boundary of the polygon.</li>
</ul>
</div>
<div id="coordinate-system" class="section level2">
<h2>Coordinate system</h2>
<p>To correctly plot spatial data, you need to know exactly what the numeric positions mean, i.e. what are they in reference to? This is called the <strong>coordinate reference system</strong> or CRS. Often spatial data is described in terms of latitude and longitude. You can check this with <code>st_is_longlat()</code>:</p>
<pre class="r"><code>st_is_longlat(nc)</code></pre>
<pre><code>## [1] TRUE</code></pre>
<p>You might think that if you know the latitude and longitude of a point, you know exactly where it is on the Earth. However, things are not quite so simple, because latitude and longitude are based on the assumption that the Earth is a smooth ellipsoid, which is not true. Because different approximations work better in differently places, most countries have their own approximation: this is called the <strong>geodetic datum</strong>, or just <strong>datum</strong> for short.</p>
<p>Take two minutes and watch this <a href="https://www.youtube.com/watch?v=xKGlMp__jog">simple explanation of the datum</a>: <iframe width="560" height="315" src="https://www.youtube.com/embed/xKGlMp__jog" frameborder="0" allow="autoplay; encrypted-media" allowfullscreen></iframe></p>
<p>To get the datum and other coordinate system metadata, use <code>st_crs()</code>:</p>
<pre class="r"><code>st_crs(nc)</code></pre>
<pre><code>## Coordinate Reference System:
## EPSG: 4267
## proj4string: "+proj=longlat +datum=NAD27 +no_defs"</code></pre>
<p>Here the datum is “NAD27”, the <a href="https://en.wikipedia.org/wiki/North_American_Datum">North American Datum</a> of 1927 (NAD27)</p>
<p>In this class, you won’t have to worry too much about the datum as sf and ggplot2 will take care of the details for you. But it’s good to know why it exists and how to identify it if something goes wrong.</p>
</div>
</div>
<div id="spatial-visualization" class="section level1">
<h1>Spatial Visualization</h1>
<div id="setup" class="section level2">
<h2>Setup</h2>
<p>Since sf is so new, support for it in ggplot2 is also very new. That means you’ll need to a recent version of ggplot2. Check that the install has succeeded by loading the tidyverse (you should have done this above) and then running <code>?geom_sf</code>. If you can’t find the documentation for <code>geom_sf</code>, something has gone wrong. The first thing to try is restarting R so that you have a clean session. Installing ggplot2 is tricky if you already have it loaded!</p>
<p>Let’s load two datasets:</p>
<pre class="r"><code>nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
states <- sf::st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))</code></pre>
<div id="geom_sf" class="section level4">
<h4><code>geom_sf()</code></h4>
<p>The easiest way to get started is to supply an sf object to <code>geom_sf()</code>:</p>
<pre class="r"><code>ggplot() +
geom_sf(data = nc)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-16-1.png" /><!-- --></p>
<p>You can also write this as:</p>
<pre class="r"><code>ggplot(nc) +
geom_sf()</code></pre>
<p>Notice that ggplot2 takes care of setting the aspect ratio correctly.</p>
<p>You can supply other aesthetics: for polygons, <code>fill</code> is most useful:</p>
<pre class="r"><code>ggplot() +
geom_sf(aes(fill = AREA), data = nc, colour = "white")</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-18-1.png" /><!-- --></p>
<p>When you include multiple layers, ggplot2 will take care of ensuring that they all have a common CRS so that it makes sense to overlay them.</p>
<pre class="r"><code>ggplot() +
geom_sf(data = states) +
geom_sf(data = nc)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-19-1.png" /><!-- --></p>
<p>You can combine <code>geom_sf()</code> with other geoms. In this case, <code>x</code> and <code>y</code> positions are assumed be in the same CRS as the sf object (typically these will be longitude and latitude).</p>
<pre class="r"><code>ggplot() +
geom_sf(data = nc) +
annotate("point", x = -80, y = 35, colour = "red", size = 4)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-20-1.png" /><!-- --></p>
</div>
<div id="coord_sf" class="section level4">
<h4><code>coord_sf()</code></h4>
<p>You’ll need to use <code>coord_sf()</code> for two reasons:</p>
<ol style="list-style-type: decimal">
<li>You want to zoom into a specified region of the plot by using <code>xlim</code> and <code>ylim</code></li>
</ol>
<pre class="r"><code>ggplot() +
geom_sf(data = nc) +
annotate("point", x = -80, y = 35, colour = "red", size = 4) +
coord_sf(xlim = c(-81, -79), ylim = c(34, 36))</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-21-1.png" /><!-- --></p>
<ol start="2" style="list-style-type: decimal">
<li>You want to override to use a specific projection. If you don’t specify the <code>crs</code> argument, it just uses the one provided in the first layer. The following example uses “USA_Contiguous_Albers_Equal_Area_Conic”. The easiest way to supply the CRS is as a EPSG ID. I found this ID (102003) with a little googling.</li>
</ol>
<pre class="r"><code>ggplot() +
geom_sf(data = states) +
coord_sf(crs = st_crs(102003))</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-22-1.png" /><!-- --></p>
</div>
</div>
<div id="reading-and-writing" class="section level2">
<h2>Reading and writing</h2>
<p>As we’ve seen above, reading spatial data from an external file can be done by</p>
<pre class="r"><code>filename <- system.file("shape/nc.shp", package="sf")
nc <- st_read(filename)</code></pre>
<pre><code>## Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/3.5/Resources/library/sf/shape/nc.shp' using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID): 4267
## proj4string: +proj=longlat +datum=NAD27 +no_defs</code></pre>
<p>we can suppress the output by adding argument <code>quiet=TRUE</code> or by using the otherwise nearly identical but more quiet</p>
<pre class="r"><code>nc <- read_sf(filename)</code></pre>
<p>Writing takes place in the same fashion, using <code>st_write</code>:</p>
<pre class="r"><code>st_write(nc, "nc.shp")</code></pre>
<p>If we repeat this, we get an error message that the file already exists, and we can overwrite by</p>
<pre class="r"><code>st_write(nc, "nc.shp", delete_layer = TRUE)</code></pre>
<p>or its quiet alternative that does this by default,</p>
<pre class="r"><code>write_sf(nc, "nc.shp") # silently overwrites</code></pre>
</div>
<div id="geometrycollection" class="section level2">
<h2>Geometrical operations</h2>
<p>The standard for simple feature access defines a number of geometrical operations.</p>
<p><code>st_is_valid</code> and <code>st_is_simple</code> return a boolean indicating whether a geometry is valid or simple.</p>
<pre class="r"><code>st_is_valid(nc[1:2,])</code></pre>
<pre><code>## [1] TRUE TRUE</code></pre>
<p><code>st_distance</code> returns a dense numeric matrix with distances between geometries. <code>st_relate</code> returns a character matrix with the <a href="https://en.wikipedia.org/wiki/DE-9IM#Illustration">DE9-IM</a> values for each pair of geometries:</p>
<pre class="r"><code>nc_nad83 = st_transform(nc, 32119) # reproject to NAD83 / North Carolina
st_distance(nc_nad83[c(1,4,22),], nc_nad83[c(1, 33,55,56),])</code></pre>
<pre><code>## Units: m
## [,1] [,2] [,3] [,4]
## [1,] 0.00 312176.2 128338.51 475608.8
## [2,] 440548.35 114938.1 590417.79 0.0
## [3,] 18943.74 352708.6 78754.75 517511.6</code></pre>
<pre class="r"><code>st_relate(nc_nad83[1:5,], nc_nad83[1:4,])</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] "2FFF1FFF2" "FF2F11212" "FF2FF1212" "FF2FF1212"
## [2,] "FF2F11212" "2FFF1FFF2" "FF2F11212" "FF2FF1212"
## [3,] "FF2FF1212" "FF2F11212" "2FFF1FFF2" "FF2FF1212"
## [4,] "FF2FF1212" "FF2FF1212" "FF2FF1212" "2FFF1FFF2"
## [5,] "FF2FF1212" "FF2FF1212" "FF2FF1212" "FF2FF1212"</code></pre>
<p>The commands <code>st_intersects</code>, <code>st_disjoint</code>, <code>st_touches</code>, <code>st_crosses</code>, <code>st_within</code>, <code>st_contains</code>, <code>st_overlaps</code>, <code>st_equals</code>, <code>st_covers</code>, <code>st_covered_by</code>, <code>st_equals_exact</code> and <code>st_is_within_distance</code> return a sparse matrix with matching (TRUE) indexes, or a full logical matrix:</p>
<pre class="r"><code>st_intersects(nc_nad83[1:5,], nc_nad83[1:4,])</code></pre>
<pre><code>## Sparse geometry binary predicate list of length 5, where the predicate was `intersects'
## 1: 1, 2
## 2: 1, 2, 3
## 3: 2, 3
## 4: 4
## 5: (empty)</code></pre>
<pre class="r"><code>st_intersects(nc_nad83[1:5,], nc_nad83[1:4,], sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE FALSE FALSE
## [2,] TRUE TRUE TRUE FALSE
## [3,] FALSE TRUE TRUE FALSE
## [4,] FALSE FALSE FALSE TRUE
## [5,] FALSE FALSE FALSE FALSE</code></pre>
<p>The commands <code>st_buffer</code>, <code>st_boundary</code>, <code>st_convexhull</code>, <code>st_union_cascaded</code>, <code>st_simplify</code>, <code>st_triangulate</code>, <code>st_polygonize</code>, <code>st_centroid</code>, <code>st_segmentize</code>, and <code>st_union</code> return new geometries, e.g.:</p>
<pre class="r"><code>sel <- c(1,5,14)
geom = st_geometry(nc_nad83)
buf <- st_buffer(geom, dist = 30000)
plot(buf, border = 'red')
plot(geom, add = TRUE)
plot(st_buffer(geom, -5000), add = TRUE, border = 'blue')</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-31-1.png" /><!-- --></p>
<p>Commands <code>st_intersection</code>, <code>st_union</code>, <code>st_difference</code>, <code>st_sym_difference</code> return new geometries that are a function of pairs of geometries:</p>
<pre class="r"><code>par(mar = rep(0,4))
u <- st_union(nc_nad83)
plot(u)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-32-1.png" /><!-- --></p>
<p>The following code shows how computing an intersection between two polygons may yield a <code>GEOMETRYCOLLECTION</code> with a point, line and polygon:</p>
<pre class="r"><code>opar <- par(mfrow = c(1, 2))
a <- st_polygon(list(cbind(c(0,0,7.5,7.5,0),c(0,-1,-1,0,0))))
b <- st_polygon(list(cbind(c(0,1,2,3,4,5,6,7,7,0),c(1,0,.5,0,0,0.5,-0.5,-0.5,1,1))))
plot(a, ylim = c(-1,1))
title("intersecting two polygons:")
plot(b, add = TRUE, border = 'red')
(i <- st_intersection(a,b))</code></pre>
<pre><code>## GEOMETRYCOLLECTION (POINT (1 0), LINESTRING (4 0, 3 0), POLYGON ((5.5 0, 7 0, 7 -0.5, 6 -0.5, 5.5 0)))</code></pre>
<pre class="r"><code>plot(a, ylim = c(-1,1))
title("GEOMETRYCOLLECTION")
plot(b, add = TRUE, border = 'red')
plot(i, add = TRUE, col = 'green', lwd = 2)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-33-1.png" /><!-- --></p>
<pre class="r"><code>par(opar)</code></pre>
</div>
<div id="non-valid-geometries" class="section level2">
<h2>Non-valid geometries</h2>
<p>Invalid geometries are for instance self-intersecting lines (left) or polygons with slivers (middle) or self-intersections (right).</p>
<pre class="r"><code>x1 <- st_linestring(cbind(c(0,1,0,1),c(0,1,1,0)))
x2 <- st_polygon(list(cbind(c(0,1,1,1,0,0),c(0,0,1,0.6,1,0))))
x3 <- st_polygon(list(cbind(c(0,1,0,1,0),c(0,1,1,0,0))))
st_is_simple(st_sfc(x1))</code></pre>
<pre><code>## [1] FALSE</code></pre>
<pre class="r"><code>st_is_valid(st_sfc(x2,x3))</code></pre>
<pre><code>## [1] FALSE FALSE</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-35-1.png" /><!-- --></p>
</div>
</div>
<div id="units" class="section level1">
<h1>Units</h1>
<p>Where possible geometric operations such as <code>st_distance()</code>, <code>st_length()</code> and <code>st_area()</code> report results with a units attribute appropriate for the CRS:</p>
<pre class="r"><code>a <- st_area(nc[1,])
attributes(a)</code></pre>
<pre><code>## $units
## $numerator
## [1] "m" "m"
##
## $denominator
## character(0)
##
## attr(,"class")
## [1] "symbolic_units"
##
## $class
## [1] "units"</code></pre>
<p>The <strong>units</strong> package can be used to convert between units:</p>
<pre class="r"><code>units::set_units(a, km^2) # result in square kilometers</code></pre>
<pre><code>## 1137.389 km^2</code></pre>
<pre class="r"><code>units::set_units(a, ha) # result in hectares</code></pre>
<pre><code>## 113738.9 ha</code></pre>
<p>The result can be stripped of their attributes if needs be:</p>
<pre class="r"><code>as.numeric(a)</code></pre>
<pre><code>## [1] 1137388604</code></pre>
<div id="geometrical-operations" class="section level2">
<h2>Geometrical operations</h2>
<p>All geometrical operations <code>st_op(x)</code> or or <code>st_op2(x,y)</code> work both for <code>sf</code> objects as well as <code>sfc</code> objects <code>x</code> and <code>y</code>; since the operations work on the geometries, the non-geometries parts of an <code>sf</code> object are simply discarded. Also, all binary operations <code>st_op2(x,y)</code> called with a single argument, as <code>st_op2(x)</code>, are handled as <code>st_op2(x,x)</code>.</p>
<p>We will illustrate the geometrical operations on a very simple dataset:</p>
<pre class="r"><code>b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
b1 = b0 + 2
b2 = b0 + c(-0.2, 2)
x = st_sfc(b0, b1, b2)
a0 = b0 * 0.8
a1 = a0 * 0.5 + c(2, 0.7)
a2 = a0 + 1
a3 = b0 * 0.5 + c(2, -0.5)
y = st_sfc(a0,a1,a2,a3)
plot(x, border = 'red')
plot(y, border = 'green', add = TRUE)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-39-1.png" /><!-- --></p>
<div id="unary-operations" class="section level3">
<h3>Unary operations</h3>
<p><code>st_is_valid</code> returns whether polygon geometries are topologically valid:</p>
<pre class="r"><code>b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
b1 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(0,-1), c(-1,-1))))
st_is_valid(st_sfc(b0,b1))</code></pre>
<pre><code>## [1] TRUE FALSE</code></pre>
<p>and <code>st_is_simple</code> whether line geometries are simple:</p>
<pre class="r"><code>s = st_sfc(st_linestring(rbind(c(0,0), c(1,1))),
st_linestring(rbind(c(0,0), c(1,1),c(0,1),c(1,0))))
st_is_simple(s)</code></pre>
<pre><code>## [1] TRUE FALSE</code></pre>
<p><code>st_area</code> returns the area of polygon geometries, <code>st_length</code> the length of line geometries:</p>
<pre class="r"><code>st_area(x)</code></pre>
<pre><code>## [1] 4 4 4</code></pre>
<pre class="r"><code>st_area(st_sfc(st_point(c(0,0))))</code></pre>
<pre><code>## [1] 0</code></pre>
<pre class="r"><code>st_length(st_sfc(st_linestring(rbind(c(0,0),c(1,1),c(1,2))), st_linestring(rbind(c(0,0),c(1,0)))))</code></pre>
<pre><code>## [1] 2.414214 1.000000</code></pre>
<pre class="r"><code>st_length(st_sfc(st_multilinestring(list(rbind(c(0,0),c(1,1),c(1,2))),rbind(c(0,0),c(1,0))))) # ignores 2nd part!</code></pre>
<pre><code>## [1] 2.414214</code></pre>
</div>
<div id="binary-operations-distance-and-relate" class="section level3">
<h3>Binary operations: distance and relate</h3>
<p><code>st_distance</code> computes the shortest distance matrix between geometries; this is a dense matrix:</p>
<pre class="r"><code>st_distance(x,y)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] 0.0000000 0.6 0 0.500000
## [2,] 0.2828427 0.0 0 1.000000
## [3,] 0.2000000 0.8 0 1.220656</code></pre>
<p><code>st_relate</code> returns a dense character matrix with the DE9-IM relationships between each pair of geometries:</p>
<pre class="r"><code>st_relate(x,y)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] "212FF1FF2" "FF2FF1212" "212101212" "FF2FF1212"
## [2,] "FF2FF1212" "212101212" "212101212" "FF2FF1212"
## [3,] "FF2FF1212" "FF2FF1212" "212101212" "FF2FF1212"</code></pre>
<p>element [i,j] of this matrix has nine characters, refering to relationship between x[i] and y[j], encoded as <span class="math inline"><em>I</em><sub><em>x</em></sub><em>I</em><sub><em>y</em></sub>, <em>I</em><sub><em>x</em></sub><em>B</em><sub><em>y</em></sub>, <em>I</em><sub><em>x</em></sub><em>E</em><sub><em>y</em></sub>, <em>B</em><sub><em>x</em></sub><em>I</em><sub><em>y</em></sub>, <em>B</em><sub><em>x</em></sub><em>B</em><sub><em>y</em></sub>, <em>B</em><sub><em>x</em></sub><em>E</em><sub><em>y</em></sub>, <em>E</em><sub><em>x</em></sub><em>I</em><sub><em>y</em></sub>, <em>E</em><sub><em>x</em></sub><em>B</em><sub><em>y</em></sub>, <em>E</em><sub><em>x</em></sub><em>E</em><sub><em>y</em></sub></span> where <span class="math inline"><em>I</em></span> refers to interior, <span class="math inline"><em>B</em></span> to boundary, and <span class="math inline"><em>E</em></span> to exterior, and e.g. <span class="math inline"><em>B</em><sub><em>x</em></sub><em>I</em><sub><em>y</em></sub></span> the dimensionality of the intersection of the the boundary <span class="math inline"><em>B</em><sub><em>x</em></sub></span> of x[i] and the interior <span class="math inline"><em>I</em><sub><em>y</em></sub></span> of y[j], which is one of {0,1,2,F}, indicating zero-, one-, two-dimension intersection, and (F) no intersection, respectively.</p>
<p><img src="04_assets/DE9-IM.png" alt="DE9-IM" /> Reading from left-to-right and top-to-bottom, the DE-9IM(a,b) string code is ‘212101212’, the compact representation of <span class="math inline"><em>I</em><sub><em>x</em></sub><em>I</em><sub><em>y</em></sub> = 2, <em>I</em><sub><em>x</em></sub><em>B</em><sub><em>y</em></sub> = 1, <em>I</em><sub><em>x</em></sub><em>E</em><sub><em>y</em></sub> = 2, <em>B</em><sub><em>x</em></sub><em>I</em><sub><em>y</em></sub> = 1, <em>B</em><sub><em>x</em></sub><em>B</em><sub><em>y</em></sub> = 0, <em>B</em><sub><em>x</em></sub><em>E</em><sub><em>y</em></sub> = 1, <em>E</em><sub><em>x</em></sub><em>I</em><sub><em>y</em></sub> = 2, <em>E</em><sub><em>x</em></sub><em>B</em><sub><em>y</em></sub> = 1, <em>E</em><sub><em>x</em></sub><em>E</em><sub><em>y</em></sub> = 2</span>. Figure from <a href="https://en.wikipedia.org/wiki/DE-9IM#Illustration">here</a>.</p>
</div>
<div id="binary-logical-operations" class="section level3">
<h3>Binary logical operations:</h3>
<p>Binary logical operations return either a sparse matrix</p>
<pre class="r"><code>st_intersects(x,y)</code></pre>
<pre><code>## Sparse geometry binary predicate list of length 3, where the predicate was `intersects'
## 1: 1, 3
## 2: 2, 3
## 3: 3</code></pre>
<p>or a dense matrix</p>
<pre class="r"><code>st_intersects(x, x, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE FALSE
## [3,] TRUE FALSE TRUE</code></pre>
<pre class="r"><code>st_intersects(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] TRUE FALSE TRUE FALSE
## [2,] FALSE TRUE TRUE FALSE
## [3,] FALSE FALSE TRUE FALSE</code></pre>
<p>where list element <code>i</code> of a sparse matrix contains the indices of the <code>TRUE</code> elements in row <code>i</code> of the the dense matrix. For large geometry sets, dense matrices take up a lot of memory and are mostly filled with <code>FALSE</code> values, hence the default is to return a sparse matrix.</p>
<p><code>st_intersects</code> returns for every geometry pair whether they intersect (dense matrix), or which elements intersect (sparse).</p>
<p>Other binary predicates include (using sparse for readability):</p>
<pre class="r"><code>st_disjoint(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] FALSE TRUE FALSE TRUE
## [2,] TRUE FALSE FALSE TRUE
## [3,] TRUE TRUE FALSE TRUE</code></pre>
<pre class="r"><code>st_touches(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE</code></pre>
<pre class="r"><code>st_crosses(s, s, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2]
## [1,] FALSE FALSE
## [2,] FALSE FALSE</code></pre>
<pre class="r"><code>st_within(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE</code></pre>
<pre class="r"><code>st_contains(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] TRUE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE</code></pre>
<pre class="r"><code>st_overlaps(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] FALSE FALSE TRUE FALSE
## [2,] FALSE TRUE TRUE FALSE
## [3,] FALSE FALSE TRUE FALSE</code></pre>
<pre class="r"><code>st_equals(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE</code></pre>
<pre class="r"><code>st_covers(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] TRUE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE</code></pre>
<pre class="r"><code>st_covered_by(x, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE</code></pre>
<pre class="r"><code>st_covered_by(y, y, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] TRUE FALSE FALSE FALSE
## [2,] FALSE TRUE FALSE FALSE
## [3,] FALSE FALSE TRUE FALSE
## [4,] FALSE FALSE FALSE TRUE</code></pre>
<pre class="r"><code>st_equals_exact(x, y,0.001, sparse = FALSE)</code></pre>
<pre><code>## [,1] [,2] [,3] [,4]
## [1,] FALSE FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE FALSE</code></pre>
</div>
<div id="operations-returning-a-geometry" class="section level3">
<h3>Operations returning a geometry</h3>
<pre class="r"><code>u = st_union(x)
plot(u)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-48-1.png" /><!-- --></p>
<pre class="r"><code>par(mfrow=c(1,2), mar = rep(0,4))
plot(st_buffer(u, 0.2))
plot(u, border = 'red', add = TRUE)
plot(st_buffer(u, 0.2), border = 'grey')
plot(u, border = 'red', add = TRUE)
plot(st_buffer(u, -0.2), add = TRUE)</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-49-1.png" /><!-- --></p>
<pre class="r"><code>plot(st_boundary(x))</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-50-1.png" /><!-- --></p>
<pre class="r"><code>par(mfrow = c(1:2))
plot(st_convex_hull(x))
plot(st_convex_hull(u))</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-51-1.png" /><!-- --></p>
<pre class="r"><code>par(mfrow = c(1,1))</code></pre>
<pre class="r"><code>par(mfrow=c(1,2))
plot(x)
plot(st_centroid(x), add = TRUE, col = 'red')
plot(x)
plot(st_centroid(u), add = TRUE, col = 'red')</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-52-1.png" /><!-- --></p>
<p>The intersection of two geometries is the geometry covered by both; it is obtained by <code>st_intersection</code>:</p>
<pre class="r"><code>plot(x)
plot(y, add = TRUE)
plot(st_intersection(st_union(x),st_union(y)), add = TRUE, col = 'red')</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-53-1.png" /><!-- --></p>
<p>To get <em>everything but</em> the intersection, use <code>st_difference</code> or st_sym_difference`:</p>
<pre class="r"><code>par(mfrow=c(2,2), mar = c(0,0,1,0))
plot(x, col = '#ff333388');
plot(y, add=TRUE, col='#33ff3388')
title("x: red, y: green")
plot(x, border = 'grey')
plot(st_difference(st_union(x),st_union(y)), col = 'lightblue', add = TRUE)
title("difference(x,y)")
plot(x, border = 'grey')
plot(st_difference(st_union(y),st_union(x)), col = 'lightblue', add = TRUE)
title("difference(y,x)")
plot(x, border = 'grey')
plot(st_sym_difference(st_union(y),st_union(x)), col = 'lightblue', add = TRUE)
title("sym_difference(x,y)")</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-54-1.png" /><!-- --></p>
<p>Function <code>st_segmentize</code> adds points to straight line sections of a lines or polygon object:</p>
<pre class="r"><code>par(mfrow=c(1,3),mar=c(1,1,0,0))
pts = rbind(c(0,0),c(1,0),c(2,1),c(3,1))
ls = st_linestring(pts)
plot(ls)
points(pts)
ls.seg = st_segmentize(ls, 0.3)
plot(ls.seg)
pts = ls.seg
points(pts)
pol = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0))))
pol.seg = st_segmentize(pol, 0.3)
plot(pol.seg, col = 'grey')
points(pol.seg[[1]])</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-55-1.png" /><!-- --></p>
<p>Function <code>st_polygonize</code> polygonizes a multilinestring, as far as the points form a closed polygon:</p>
<pre class="r"><code>par(mfrow=c(1,2),mar=c(0,0,1,0))
mls = st_multilinestring(list(matrix(c(0,0,0,1,1,1,0,0),,2,byrow=TRUE)))
x = st_polygonize(mls)
plot(mls, col = 'grey')
title("multilinestring")
plot(x, col = 'grey')
title("polygon")</code></pre>
<p><img src="04_Spatial_with_sf_files/figure-html/unnamed-chunk-56-1.png" /><!-- --></p>
<p>Further reading: 1. Much of the material below was taken from the <a href="https://cran.r-project.org/web/packages/sf/vignettes">sf vignettes available here</a>. 2. S. Scheider, B. Gräler, E. Pebesma, C. Stasch, 2016. Modelling spatio-temporal information generation. Int J of Geographic Information Science, 30 (10), 1980-2008. (<a href="http://pebesma.staff.ifgi.de/generativealgebra.pdf">pdf</a>) 3. Stasch, C., S. Scheider, E. Pebesma, W. Kuhn, 2014. Meaningful Spatial Prediction and Aggregation. Environmental Modelling & Software, 51, (149–165, <a href="http://dx.doi.org/10.1016/j.envsoft.2013.09.006">open access</a>).</p>
</div>
</div>
</div>