-
Notifications
You must be signed in to change notification settings - Fork 303
/
Copy pathplot.R
1070 lines (997 loc) · 38.1 KB
/
plot.R
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
kw_dflt = function(x, key.pos) {
if (is.null(key.pos) || key.pos[1] == 0) # no key:
return(lcm(0))
font_scale = par("ps") / 12
if (key.pos[1] == -1)
lcm(1.8 * font_scale)
else if (key.pos[1] %in% c(2, 4) && (is.character(x[[1]]) || is.factor(x[[1]]))) {
strings = if (is.factor(x[[1]]))
levels(x[[1]])
else
x[[1]]
lcm(cm(max(strwidth(strings, "inches"))) * 1.3 + font_scale) # cm
#lcm(cm(max(strwidth(strings, "inches"))) * 1.3) # cm
} else
lcm(1.8 * font_scale)
}
# like cut.default, but only return integers, and allow for duplicate breaks:
sf_cut = function(values, breaks, include.lowest = TRUE) {
i = findInterval(values, breaks, left.open = TRUE)
if (include.lowest)
i[values == breaks[1]] = 1
i[i == 0 | i == length(breaks)] = NA
i
}
#' plot sf object
#'
#' plot one or more attributes of an sf object on a map
#' Plot sf object
#'
#' @param x object of class sf
#' @param y ignored
#' @param ... further specifications, see \link{plot_sf} and \link{plot} and details.
#' @param main title for plot (\code{NULL} to remove)
#' @param pal palette function, similar to \link{rainbow}, or palette values; if omitted, \code{sf.colors} is used
#' @param nbreaks number of colors breaks (ignored for \code{factor} or \code{character} variables)
#' @param breaks either a numeric vector with the actual breaks, or a name of a method accepted by the \code{style} argument of \link[classInt]{classIntervals}
#' @param max.plot integer; lower boundary to maximum number of attributes to plot; the default value (9) can be overridden by setting the global option \code{sf_max.plot}, e.g. \code{options(sf_max.plot=2)}
#' @param key.pos numeric; side to plot a color key: 1 bottom, 2 left, 3 top, 4 right; set to \code{NULL} to omit key completely, 0 to only not plot the key, or -1 to select automatically. If multiple columns are plotted in a single function call by default no key is plotted and every submap is stretched individually; if a key is requested (and \code{col} is missing) all maps are colored according to a single key. Auto select depends on plot size, map aspect, and, if set, parameter \code{asp}. If it has lenght 2, the second value, ranging from 0 to 1, determines where the key is placed in the available space (default: 0.5, center).
#' @param key.width amount of space reserved for the key (incl. labels), thickness/width of the scale bar
#' @param key.length amount of space reserved for the key along its axis, length of the scale bar
#' @param pch plotting symbol
#' @param cex symbol size
#' @param bg symbol background color
#' @param lty line type
#' @param lwd line width
#' @param col color for plotting features; if \code{length(col)} does not equal 1 or \code{nrow(x)}, a warning is emitted that colors will be recycled. Specifying \code{col} suppresses plotting the legend key.
#' @param border color of polygon border(s); using \code{NA} hides them
#' @param add logical; add to current plot? Note that when using \code{add=TRUE}, you may have to set \code{reset=FALSE} in the first plot command.
#' @param type plot type: 'p' for points, 'l' for lines, 'b' for both
#' @param reset logical; if \code{FALSE}, keep the plot in a mode that allows adding further map elements; if \code{TRUE} restore original mode after plotting \code{sf} objects with attributes; see details.
#' @param logz logical; if \code{TRUE}, use log10-scale for the attribute variable. In that case, \code{breaks} and \code{at} need to be given as log10-values; see examples.
#' @param extent object with an \code{st_bbox} method to define plot extent; defaults to \code{x}
#' @param xlim numeric; x-axis limits; overrides \code{extent}
#' @param ylim numeric; y-axis limits; overrides \code{extent}
#' @param compact logical; compact sub-plots over plotting space?
#' @method plot sf
#' @name plot
#' @details \code{plot.sf} maximally plots \code{max.plot} maps with colors following from attribute columns,
#' one map per attribute. It uses \code{sf.colors} for default colors. For more control over placement of individual maps,
#' set parameter \code{mfrow} with \link{par} prior to plotting, and plot single maps one by one; note that this only works
#' in combination with setting parameters \code{key.pos=NULL} (no legend) and \code{reset=FALSE}.
#'
#' \code{plot.sfc} plots the geometry, additional parameters can be passed on
#' to control color, lines or symbols.
#'
#' When setting \code{reset} to \code{FALSE}, the original device parameters are lost, and the device must be reset using \code{dev.off()} in order to reset it.
#'
#' parameter \code{at} can be set to specify where labels are placed along the key; see examples.
#'
#' The features are plotted in the order as they apppear in the sf object. See examples for when a different plotting order is wanted.
#'
#' @examples
#' nc = st_read(system.file("gpkg/nc.gpkg", package="sf"), quiet = TRUE)
#' # plot single attribute, auto-legend:
#' plot(nc["SID74"])
#' # plot multiple:
#' plot(nc[c("SID74", "SID79")]) # better use ggplot2::geom_sf to facet and get a single legend!
#' # adding to a plot of an sf object only works when using reset=FALSE in the first plot:
#' plot(nc["SID74"], reset = FALSE)
#' plot(st_centroid(st_geometry(nc)), add = TRUE)
#' # log10 z-scale:
#' plot(nc["SID74"], logz = TRUE, breaks = c(0,.5,1,1.5,2), at = c(0,.5,1,1.5,2))
#' # and we need to reset the plotting device after that, e.g. by
#' layout(1)
#' # when plotting only geometries, the reset=FALSE is not needed:
#' plot(st_geometry(nc))
#' plot(st_geometry(nc)[1], col = 'red', add = TRUE)
#' # add a custom legend to an arbitray plot:
#' layout(matrix(1:2, ncol = 2), widths = c(1, lcm(2)))
#' plot(1)
#' .image_scale(1:10, col = sf.colors(9), key.length = lcm(8), key.pos = 4, at = 1:10)
#' # manipulate plotting order, plot largest polygons first:
#' p = st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
#' x = st_sf(a=1:4, st_sfc(p, p * 2, p * 3, p * 4)) # plot(x, col=2:5) only shows the largest polygon!
#' plot(x[order(st_area(x), decreasing = TRUE),], col = 2:5) # plot largest polygons first
#'
#' @export
plot.sf <- function(x, y, ..., main, pal = NULL, nbreaks = 10, breaks = "pretty",
max.plot = getOption("sf_max.plot", default = 9),
key.pos = get_key_pos(x, ...), key.length = .618, key.width = kw_dflt(x, key.pos),
reset = TRUE, logz = FALSE, extent = x, xlim = st_bbox(extent)[c(1,3)],
ylim = st_bbox(extent)[c(2,4)], compact = FALSE) {
stopifnot(missing(y))
nbreaks.missing = missing(nbreaks)
key.pos.missing = missing(key.pos)
max_plot_missing = missing(max.plot)
dots = list(...)
col_missing = is.null(dots$col)
breaks_numeric = is.numeric(breaks)
reset_layout_needed = reset
x = swap_axes_if_needed(x)
opar = par(no.readonly = TRUE)
if (ncol(x) > 2 && !isTRUE(dots$add)) { # multiple maps to plot...
cols = setdiff(names(x), attr(x, "sf_column"))
lt = .get_layout(st_bbox(x), min(max.plot, length(cols)), par("din"), key.pos[1], key.width)
if (key.pos.missing || key.pos[1] == -1)
key.pos = lt$key.pos
layout(lt$m, widths = lt$widths, heights = lt$heights, respect = compact)
if (isTRUE(dots$axes))
par(mar = c(2.1, 2.1, 1.2, 0))
else
par(mar = c(0, 0, 1.2, 0))
if (max_plot_missing)
max.plot = prod(lt$mfrow)
if (isTRUE(is.finite(max.plot)) && ncol(x) - 1 > max.plot &&
max_plot_missing && is.null(options("sf_max.plot")[[1]]))
warning(paste("plotting the first", max.plot, "out of", ncol(x)-1,
"attributes; use max.plot =", ncol(x) - 1, "to plot all"), call. = FALSE)
# col selection may have changed; set cols again:
cols = setdiff(names(x), attr(x, "sf_column"))
if (length(cols) > max.plot)
cols = cols[1:max.plot]
if (!is.null(key.pos)) {
values = do.call(c, as.data.frame(x)[cols])
if (is.character(values))
values = as.factor(values)
if (logz)
values = log10(values)
if (is.character(breaks) && is.numeric(values)) { # compute breaks from values:
v0 = values[!is.na(values)]
n.unq = length(unique(v0))
breaks = if (! all(is.na(values)) && n.unq > 1)
classInt::classIntervals(v0, min(nbreaks, n.unq),
breaks, warnSmallN = FALSE)$brks
else
range(values, na.rm = TRUE) # lowest and highest!
nbreaks = length(breaks) - 1
}
}
if (nbreaks.missing && is.numeric(breaks))
nbreaks = length(breaks) - 1
# loop over each map to plot:
lapply(cols, function(cname) plot(x[, cname], main = cname,
pal = pal, nbreaks = nbreaks, breaks = breaks, key.pos = NULL, reset = FALSE,
logz = logz, xlim = xlim, ylim = ylim,...))
for (i in seq_len(prod(lt$mfrow) - length(cols))) # empty panels:
plot.new()
# plot key?
if (!is.null(key.pos) && key.pos[1] != 0 && col_missing) {
if (is.null(pal))
pal = function(n) sf.colors(n, categorical = is.factor(values))
colors = if (is.function(pal))
pal(nbreaks)
else
pal
if (is.factor(values))
.image_scale_factor(levels(values), colors, key.pos = key.pos,
key.width = key.width, key.length = key.length, ...)
else
.image_scale(values, colors, breaks = breaks, key.pos = key.pos,
key.length = key.length, logz = logz, ...)
}
} else { # single map, or dots$add == TRUE:
if (ncol(x) == 1) { # no attributes to choose colors from: plot geometry
plot(st_geometry(x), xlim = xlim, ylim = ylim, ...)
reset_layout_needed = FALSE
} else { # generate plot with colors and possibly key
if (ncol(x) > 2) { # add = TRUE
warning("ignoring all but the first attribute")
x = x[,1]
}
# store attribute in "values":
values = x[[setdiff(names(x), attr(x, "sf_column"))]]
if (is.list(values))
stop("plotting list-columns not supported") # nocov
if (is.character(values))
values = as.factor(values)
else if (logz)
values = log10(as.numeric(values))
if (is.null(pal))
pal = function(n) sf.colors(n, categorical = is.factor(values))
else if (! col_missing)
stop("specify only one of `col' and `pal'")
if (col_missing) { # compute colors from values:
col = if (is.factor(values)) {
if (key.pos.missing && nlevels(values) > 30) # doesn't make sense:
key.pos = NULL
colors = if (is.function(pal))
pal(nlevels(values))
else
pal
colors[as.numeric(values)]
} else {
if (! inherits(values, c("POSIXt", "Date")))
values = as.numeric(values) # drop units, if any
if (is.character(breaks)) { # compute breaks from values:
v0 = values[!is.na(values)]
n.unq = length(unique(v0))
breaks = if (! all(is.na(values)) && n.unq > 1)
classInt::classIntervals(v0, min(nbreaks, n.unq),
breaks, warnSmallN = FALSE)$brks
else
range(values, na.rm = TRUE) # lowest and highest!
}
# this is necessary if breaks were specified either as character or as numeric
# "pretty" takes nbreaks as advice only:
nbreaks = length(breaks) - 1
cuts = if (all(is.na(values)))
rep(NA_integer_, length(values))
else if (!breaks_numeric && diff(range(values, na.rm = TRUE)) == 0)
ifelse(is.na(values), NA_integer_, 1L)
else if (inherits(values, c("POSIXt", "Date")))
cut(values, breaks, include.lowest = TRUE)
else
sf_cut(values, breaks, include.lowest = TRUE)
colors = if (is.function(pal))
pal(nbreaks)
else
pal
colors[cuts]
}
} else {
col = dots$col
if (length(col) != 1 && length(col) != nrow(x))
warning("col is not of length 1 or nrow(x): colors will be recycled; use pal to specify a color palette")
key.pos = NULL # no key!
}
if (!isTRUE(dots$add) && !is.null(key.pos) && !all(is.na(values)) &&
(is.factor(values) || length(unique(na.omit(values))) > 1 || breaks_numeric) && # 2065
length(col) > 1) { # plot key?
switch(key.pos[1],
layout(matrix(c(2,1), nrow = 2, ncol = 1),
widths = 1, heights = c(1, key.width)), # 1 bottom
layout(matrix(c(1,2), nrow = 1, ncol = 2),
widths = c(key.width, 1), heights = 1), # 2 left
layout(matrix(c(1,2), nrow = 2, ncol = 1),
widths = 1, heights = c(key.width, 1)), # 3 top
layout(matrix(c(2,1), nrow = 1, ncol = 2),
widths = c(1, key.width), heights = 1) # 4 right
)
if (is.factor(values)) {
.image_scale_factor(levels(values), colors, key.pos = key.pos,
key.width = key.width, key.length = key.length, ...)
} else
.image_scale(values, colors, breaks = breaks, key.pos = key.pos,
key.length = key.length, logz = logz, ...)
} else
reset_layout_needed = FALSE # as we didn't call layout()
# plot the map:
if (!isTRUE(dots$add)) {
mar = c(1, 1, 1.2, 1)
if (isTRUE(dots$axes))
mar[1:2] = 2.1
par(mar = mar)
}
if (col_missing)
plot(st_geometry(x), col = col, xlim = xlim, ylim = ylim, ...)
else
plot(st_geometry(x), xlim = xlim, ylim = ylim, ...)
}
if (! isTRUE(dots$add)) { # title?
if (missing(main)) {
main = setdiff(names(x), attr(x, "sf_column"))
if (length(main) && inherits(x[[main]], "units"))
main = make_unit_label(main, x[[main]])
}
localTitle <- function(..., extent, col, bg, pch, cex, lty, lwd, axes, type, bgMap,
border, graticule, xlim, ylim, asp, bgc, xaxs, yaxs, lab, setParUsrBB,
expandBB, col_graticule, at, lon, lat, crs, datum, ndiscr, margin) # absorb
title(...)
localTitle(main, ...)
}
}
if (!isTRUE(dots$add) && reset) { # reset device:
if (reset_layout_needed)
layout(matrix(1))
par(opar)
}
invisible()
}
swap_axes_if_needed = function(x) {
crs = st_crs(x)
if (st_axis_order() && !is.na(crs) && crs$yx)
st_transform(x, pipeline = "+proj=pipeline +step +proj=axisswap +order=2,1")
else
x
}
#' @name plot
#' @export
get_key_pos = function(x, ...) {
bb = st_bbox(x)
if (anyNA(bb) || (inherits(x, "sf") && ncol(x) > 2))
NULL
else {
pin = par("pin") # (width, height)
asp_plt = pin[2]/pin[1] # y/x: < 1 means wide
asp_box = diff(bb[c(4,2)]) / diff(bb[c(3,1)])
asp = list(...)$asp
if (is.null(asp))
asp <- ifelse(isTRUE(st_is_longlat(x)), 1/cos((mean(bb[c(2,4)]) * pi)/180), 1.0)
asp_box = asp_box * asp
if (!is.finite(asp_box) || asp_box < asp_plt) # plot is wider than device: below
1
else # plot is taller than device: to the right
4
}
}
#' @name plot
#' @method plot sfc_POINT
#' @export
plot.sfc_POINT = function(x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1, lty = 1,
type = 'p', add = FALSE) {
stopifnot(missing(y))
if (! add)
plot_sf(x, ...)
npts = length(x)
pch = rep(pch, length.out = npts)
col = rep(col, length.out = npts)
bg = rep(bg, length.out = npts)
cex = rep(cex, length.out = npts)
mat = t(matrix(unlist(x, use.names = FALSE), ncol = length(x))) #933
if (!is.null(mat)) {
ne = !is.na(rowMeans(mat)) ## faster than apply; #933
points(mat[ne,, drop = FALSE], pch = pch[ne], col = col[ne], bg = bg[ne],
cex = cex[ne], lwd = lwd, lty = lty, type = type)
}
}
#' @name plot
#' @method plot sfc_MULTIPOINT
#' @export
plot.sfc_MULTIPOINT = function(x, y, ..., pch = 1, cex = 1, col = 1, bg = 0, lwd = 1, lty = 1,
type = 'p', add = FALSE) {
stopifnot(missing(y))
if (! add)
plot_sf(x, ...)
n = length(x)
pch = rep(pch, length.out = n)
col = rep(col, length.out = n)
bg = rep(bg, length.out = n)
cex = rep(cex, length.out = n)
lwd = rep(lwd, length.out = n)
lty = rep(lty, length.out = n)
non_empty = ! st_is_empty(x)
lapply(seq_along(x), function(i)
if (non_empty[i])
points(x[[i]], pch = pch[i], col = col[i], bg = bg[i],
cex = cex[i], lwd = lwd[i], lty = lty[i], type = type))
invisible(NULL)
}
#' @name plot
#' @method plot sfc_LINESTRING
#' @export
plot.sfc_LINESTRING = function(x, y, ..., lty = 1, lwd = 1, col = 1, pch = 1, type = 'l',
add = FALSE) {
# FIXME: take care of lend, ljoin, and lmitre
stopifnot(missing(y))
if (! add)
plot_sf(x, ...)
lty = rep(lty, length.out = length(x))
lwd = rep(lwd, length.out = length(x))
col = rep(col, length.out = length(x))
pch = rep(pch, length.out = length(x))
non_empty = ! st_is_empty(x)
lapply(seq_along(x), function(i)
if (non_empty[i])
lines(x[[i]], lty = lty[i], lwd = lwd[i], col = col[i], pch = pch[i], type = type))
invisible(NULL)
}
#' @name plot
#' @method plot sfc_CIRCULARSTRING
#' @export
plot.sfc_CIRCULARSTRING = function(x, y, ...) {
plot(st_cast(x, "LINESTRING"), ...)
}
#' @name plot
#' @method plot sfc_MULTILINESTRING
#' @export
plot.sfc_MULTILINESTRING = function(x, y, ..., lty = 1, lwd = 1, col = 1, pch = 1, type = 'l',
add = FALSE) {
# FIXME: take care of lend, ljoin, and lmitre
stopifnot(missing(y))
if (! add)
plot_sf(x, ...)
lty = rep(lty, length.out = length(x))
lwd = rep(lwd, length.out = length(x))
col = rep(col, length.out = length(x))
pch = rep(pch, length.out = length(x))
non_empty = ! st_is_empty(x)
lapply(seq_along(x), function(i)
if (non_empty[i])
lapply(x[[i]], function(L)
lines(L, lty = lty[i], lwd = lwd[i], col = col[i], pch = pch[i], type = type)))
invisible(NULL)
}
# sf (list) -> polypath (mtrx) : rbind polygon rings with NA rows inbetween
p_bind = function(lst) {
if (length(lst) == 1)
lst[[1]]
else {
ret = vector("list", length(lst) * 2 - 1)
ret[seq(1, length(lst) * 2 - 1, by = 2)] = lst # odd elements
ret[seq(2, length(lst) * 2 - 1, by = 2)] = NA # even elements
do.call(rbind, ret) # replicates the NA to form an NA row
}
}
#' @name plot
#' @param rule see \link[graphics]{polypath}; for \code{winding}, exterior ring direction should be opposite that of the holes; with \code{evenodd}, plotting is robust against misspecified ring directions
#' @param xpd see \link[graphics]{par}; sets polygon clipping strategy; only implemented for POLYGON and MULTIPOLYGON
#' @export
plot.sfc_POLYGON = function(x, y, ..., lty = 1, lwd = 1, col = NA, cex = 1, pch = NA, border = 1,
add = FALSE, rule = "evenodd", xpd = par("xpd")) {
# FIXME: take care of lend, ljoin, and lmitre
stopifnot(missing(y))
if (! add)
plot_sf(x, ...)
lty = rep(lty, length.out = length(x))
lwd = rep(lwd, length.out = length(x))
col = rep(col, length.out = length(x))
border = rep(border, length.out = length(x))
non_empty = ! st_is_empty(x)
lapply(seq_along(x), function(i)
if (non_empty[i])
polypath(p_bind(x[[i]]), border = border[i], lty = lty[i], lwd = lwd[i], col = col[i], rule = rule, xpd = xpd))
# if (any(!is.na(pch))) {
# pch = rep(pch, length.out = length(x))
# cex = rep(cex, length.out = length(x))
# lapply(seq_along(x), function(i)
# if (non_empty[i])
# points(p_bind(x[[i]]), pch = pch[i], cex = cex[i], type = 'p'))
# }
invisible(NULL)
}
#' @name plot
#' @method plot sfc_MULTIPOLYGON
#' @export
plot.sfc_MULTIPOLYGON = function(x, y, ..., lty = 1, lwd = 1, col = NA, border = 1, add = FALSE,
rule = "evenodd", xpd = par("xpd")) {
# FIXME: take care of lend, ljoin, and lmitre
stopifnot(missing(y))
if (! add)
plot_sf(x, ...)
lty = rep(lty, length.out = length(x))
lwd = rep(lwd, length.out = length(x))
col = rep(col, length.out = length(x))
border = rep(border, length.out = length(x))
non_empty = ! st_is_empty(x)
lapply(seq_along(x), function(i)
if (non_empty[i])
lapply(x[[i]], function(L)
polypath(p_bind(L), border = border[i], lty = lty[i], lwd = lwd[i], col = col[i], rule = rule, xpd = xpd)))
invisible(NULL)
}
# plot single geometrycollection:
plot_gc = function(x, pch, cex, bg, border = 1, lty, lwd, col, add) {
lapply(x, function(subx) {
args = list(st_sfc(subx), pch = pch, cex = cex, bg = bg, border = border,
lty = lty, lwd = lwd, col = col, add = TRUE)
fn = switch(class(subx)[2],
POINT = plot.sfc_POINT,
MULTIPOINT = plot.sfc_MULTIPOINT,
LINESTRING = plot.sfc_LINESTRING,
MULTILINESTRING = plot.sfc_MULTILINESTRING,
POLYGON = plot.sfc_POLYGON,
CIRCULARSTRING = plot.sfc_CIRCULARSTRING,
MULTIPOLYGON = plot.sfc_MULTIPOLYGON,
MULTISURFACE = plot.sfc_GEOMETRYCOLLECTION,
CURVEPOLYGON = plot.sfc_GEOMETRYCOLLECTION,
COMPOUNDCURVE = plot.sfc_GEOMETRYCOLLECTION,
GEOMETRYCOLLECTION = plot.sfc_GEOMETRYCOLLECTION,
stop(paste("plotting of", class(x)[2], "not yet supported: use st_cast?"))
)
do.call(fn, args)
})
invisible(NULL)
}
#' @name plot
#' @method plot sfc_GEOMETRYCOLLECTION
#' @export
plot.sfc_GEOMETRYCOLLECTION = function(x, y, ..., pch = 1, cex = 1, bg = 0, lty = 1, lwd = 1,
col = 1, border = 1, add = FALSE) {
# FIXME: take care of lend, ljoin, xpd, and lmitre
stopifnot(missing(y))
if (! add)
plot_sf(x, ...)
cex = rep(cex, length.out = length(x))
pch = rep(pch, length.out = length(x))
lty = rep(lty, length.out = length(x))
lwd = rep(lwd, length.out = length(x))
col = rep(col, length.out = length(x))
border = rep(border, length.out = length(x))
lapply(seq_along(x), function(i) plot_gc(x[[i]],
pch = pch[i], cex = cex[i], bg = bg[i], border = border[i], lty = lty[i],
lwd = lwd[i], col = col[i]))
invisible(NULL)
}
#' @name plot
#' @method plot sfc_GEOMETRY
#' @export
plot.sfc_GEOMETRY = function(x, y, ..., pch = 1, cex = 1, bg = 0, lty = 1, lwd = 1,
col = ifelse(st_dimension(x) == 2, NA, 1), border = 1, add = FALSE) {
stopifnot(missing(y))
if (! add)
plot_sf(x, ...)
cex = rep(cex, length.out = length(x))
pch = rep(pch, length.out = length(x))
lty = rep(lty, length.out = length(x))
lwd = rep(lwd, length.out = length(x))
col = rep(col, length.out = length(x))
border = rep(border, length.out = length(x))
lapply(seq_along(x), function(i) plot_gc(st_sfc(x[[i]]),
pch = pch[i], cex = cex[i], bg = bg[i], border = border[i], lty = lty[i],
lwd = lwd[i], col = col[i]))
invisible(NULL)
}
#' @name plot
#' @method plot sfg
#' @export
plot.sfg = function(x, ...) {
plot(st_sfc(x), ...)
}
# set up plotting area & axes; reuses sp:::plot.Spatial
#' @name plot
#' @param xlim see \link{plot.window}
#' @param ylim see \link{plot.window}
#' @param asp see below, and see \link{par}
#' @param axes logical; should axes be plotted? (default FALSE)
#' @param bgc background color
#' @param xaxs see \link{par}
#' @param yaxs see \link{par}
#' @param lab see \link{par}
#' @param setParUsrBB default FALSE; set the \code{par} \dQuote{usr} bounding box; see below
#' @param bgMap object of class \code{ggmap}, or returned by function \code{RgoogleMaps::GetMap}
#' @param expandBB numeric; fractional values to expand the bounding box with,
#' in each direction (bottom, left, top, right)
#' @param graticule logical, or object of class \code{crs} (e.g., \code{st_crs(4326)} for a WGS84 graticule), or object created by \link{st_graticule}; \code{TRUE} will give the WGS84 graticule
#' or object returned by \link{st_graticule}
#' @param col_graticule color to used for the graticule (if present)
#' @export
#' @details \code{plot_sf} sets up the plotting area, axes, graticule, or webmap background; it
#' is called by all \code{plot} methods before anything is drawn.
#'
#' The argument \code{setParUsrBB} may be used to pass the logical value \code{TRUE} to functions within \code{plot.Spatial}. When set to \code{TRUE}, par(\dQuote{usr}) will be overwritten with \code{c(xlim, ylim)}, which defaults to the bounding box of the spatial object. This is only needed in the particular context of graphic output to a specified device with given width and height, to be matched to the spatial object, when using par(\dQuote{xaxs}) and par(\dQuote{yaxs}) in addition to \code{par(mar=c(0,0,0,0))}.
#'
#' The default aspect for map plots is 1; if however data are not
#' projected (coordinates are long/lat), the aspect is by default set to
#' 1/cos(My * pi/180) with My the y coordinate of the middle of the map
#' (the mean of \code{ylim}, which defaults to the y range of bounding box). This
#' implies an \href{https://en.wikipedia.org/wiki/Equirectangular_projection}{Equirectangular projection}.
#'
plot_sf = function(x, xlim = NULL, ylim = NULL, asp = NA, axes = FALSE, bgc = par("bg"), ...,
xaxs, yaxs, lab, setParUsrBB = FALSE, bgMap = NULL, expandBB = c(0,0,0,0), graticule = NA_crs_,
col_graticule = 'grey', border, extent = x) {
# sp's bbox: matrix
# min max
# x
# y
bbox = matrix(st_bbox(extent), 2, dimnames = list(c("x", "y"), c("min", "max")))
# expandBB: 1=below, 2=left, 3=above and 4=right.
expBB = function(lim, expand) c(lim[1] - expand[1] * diff(lim), lim[2] + expand[2] * diff(lim))
if (is.null(xlim))
xlim <- expBB(bbox[1,], expandBB[c(2,4)])
if (is.null(ylim))
ylim <- expBB(bbox[2,], expandBB[c(1,3)])
if (is.na(asp))
asp <- ifelse(isTRUE(st_is_longlat(x)), 1/cos((mean(ylim) * pi)/180), 1.0)
if (anyNA(bbox))
stop("NA value(s) in bounding box. Trying to plot empty geometries?")
plot.new()
args = list(xlim = xlim, ylim = ylim, asp = asp)
if (!missing(xaxs)) args$xaxs = xaxs
if (!missing(yaxs)) args$yaxs = yaxs
if (!missing(lab)) args$lab = lab
do.call(plot.window, args)
if (setParUsrBB)
par(usr=c(xlim, ylim))
pl_reg <- par("usr")
rect(xleft = pl_reg[1], ybottom = pl_reg[3], xright = pl_reg[2],
ytop = pl_reg[4], col = bgc, border = FALSE)
linAxis = function(side, ..., lon, lat, ndiscr, reset, at) axis(side = side, ...)
if (! missing(graticule)) {
g = if (isTRUE(graticule))
st_graticule(pl_reg[c(1,3,2,4)], st_crs(x), st_crs(4326), ...)
else if (inherits(graticule, "crs") && !is.na(graticule))
st_graticule(pl_reg[c(1,3,2,4)], st_crs(x), graticule, ...)
else
graticule
plot(st_geometry(g), col = col_graticule, add = TRUE)
box()
if (axes) {
sel = g$type == "E" & g$plot12
linAxis(1L, g$x_start[sel], parse(text = g$degree_label[sel]), ...)
sel = g$type == "N" & g$plot12
linAxis(2L, g$y_start[sel], parse(text = g$degree_label[sel]), ...)
}
} else if (axes) {
box()
if (isTRUE(st_is_longlat(x))) {
local_degAxis = function(side, ..., at) .degAxis(side, ...) # absorb at
local_degAxis(1, ...)
local_degAxis(2, ...)
} else {
linAxis(1, ...)
linAxis(2, ...)
}
}
localTitle <- function(..., col, bgc, pch, cex, lty, lwd, lon, lat, ndiscr, at, labels, reset) title(...)
localTitle(...)
if (!is.null(bgMap)) {
mercator = FALSE
if (inherits(bgMap, "ggmap")) {
bb = bb2merc(bgMap, "ggmap")
mercator = TRUE
} else if (all(c("lat.center","lon.center","zoom","myTile","BBOX") %in% names(bgMap))) {
# an object returned by RgoogleMaps::GetMap
bb = bb2merc(bgMap, "RgoogleMaps")
bgMap = bgMap$myTile
mercator = TRUE
} else
bb = c(xlim[1], ylim[1], xlim[2], ylim[2]) # can be any crs!
if (mercator && st_crs(x) != st_crs(3857))
warning("crs of plotting object differs from that of bgMap, which is assumed to be st_crs(3857)") # nocov
rasterImage(bgMap, bb[1], bb[2], bb[3], bb[4], interpolate = FALSE)
}
}
#' @param n integer; number of colors
#' @param cutoff.tails numeric, in `[0,0.5]` start and end values
#' @param alpha numeric, in `[0,1]`, transparency
#' @param categorical logical; do we want colors for a categorical variable? (see details)
#' @name plot
#' @export
#' @details non-categorical colors from \code{sf.colors} were taken from \link[sp]{bpy.colors}, with modified \code{cutoff.tails} defaults
#' If categorical is \code{TRUE}, default colors are from \url{https://colorbrewer2.org/} (if n < 9, Set2, else Set3).
#' @examples
#' sf.colors(10)
sf.colors = function (n = 10, cutoff.tails = c(0.35, 0.2), alpha = 1, categorical = FALSE) {
if (categorical) {
cb = if (n <= 8)
# 8-class Set2:
c('#66c2a5','#fc8d62','#8da0cb','#e78ac3','#a6d854','#ffd92f','#e5c494','#b3b3b3')
# 12-class Set3:
else c('#8dd3c7','#ffffb3','#bebada','#fb8072','#80b1d3','#fdb462','#b3de69','#fccde5','#d9d9d9','#bc80bd','#ccebc5','#ffed6f')
# TODO: deal with alpha
if (alpha != 1.0)
cb = paste0(cb, as.hexmode(ceiling(alpha * 255)))
rep(cb, length.out = n)
} else {
i = seq(0.5 * cutoff.tails[1], 1 - 0.5 * cutoff.tails[2], length.out = n)
r = ifelse(i < .25, 0, ifelse(i < .57, i / .32 - .78125, 1))
g = ifelse(i < .42, 0, ifelse(i < .92, 2 * i - .84, 1))
b = ifelse(i < .25, 4 * i, ifelse(i < .42, 1,
ifelse(i < .92, -2 * i + 1.84, i / .08 - 11.5)))
rgb(r, g, b, alpha)
}
}
# Add text to an existing (base) plot
#
#' @param labels character, text to draw (one per row of input)
#' @name plot
#' @export
#' @details `text.sf` adds text to an existing base graphic. Text is placed at the centroid of
#' each feature in \code{x}. Provide POINT features for further control of placement.
#' `points.sf` adds point symbols to an existing base graphic. If points of text are not shown
#' correctly, try setting argument `reset` to `FALSE` in the `plot()` call.
#' @examples
#' text(nc, labels = substring(nc$NAME,1,1))
text.sf = function(x, labels = row.names(x), ...) {
text(st_geometry(x), labels = labels, ...)
}
#' @export
#' @name plot
#' @param of_largest_polygon logical, passed on to \link{st_centroid}
text.sfc = function(x, labels = seq_along(x), ..., of_largest_polygon = FALSE){
x = st_centroid(x, of_largest_polygon = of_largest_polygon)
xy = st_coordinates(x)
text(xy[,1], xy[,2], labels = labels, ...)
}
#' @name plot
#' @export
points.sf = function(x, ...) {
points(st_geometry(x), ...)
}
#' @name plot
#' @export
points.sfc = function(x, ..., of_largest_polygon = FALSE) {
x = st_centroid(x, of_largest_polygon = of_largest_polygon)
xy = st_coordinates(x)
points(xy[,1], xy[,2], ...)
}
# get the aspect ratio of a bounding box, for geodetic coords true scale at mid latitude:
get_asp = function(bb) {
asp = diff(bb[c(2,4)])/diff(bb[c(1,3)])
if (!is.finite(asp)) # 0/0
asp = 1
if (isTRUE(st_is_longlat(bb)))
asp = asp / cos(mean(bb[c(2,4)]) * pi /180)
asp
}
#' functions only exported to be used internally by stars
#' @keywords internal
#' @name stars
#' @param bb ignore
#' @param n ignore
#' @param total_size ignore
#' @param key.width ignore
#' @param key.length ignore
#' @param mfrow length-2 integer vector with number of rows, columns
#' @param main main or sub title
#' @export
.get_layout = function(bb, n, total_size, key.pos, key.width, mfrow = NULL, main = NULL) {
# return list with "m" matrix, "key.pos", "widths" and "heights" fields
# if key.pos = -1 on input, it will be a return value, "optimally" placed
asp = get_asp(bb)
strip = if (is.character(main))
# strheight(main, "inches")
par("cin")[2]
else
0.0
size = function(nrow, n, asp, strip = 0) { # given nrow n asp, what size does a single tile occupy?
ncol = ceiling(n / nrow)
xsize = total_size[1] / ncol
ysize = xsize * asp + strip
if (xsize * ysize * n > prod(total_size)) {
ysize = total_size[2] / nrow - strip
xsize = ysize / asp
}
c(xsize, ysize)
}
sz = vapply(1:n, function(nrow) size(nrow, n, asp, strip), c(0.0, 0.0))
if (is.null(mfrow)) {
nrow = which.max(apply(sz, 2, prod))
ncol = ceiling(n / nrow)
} else {
stopifnot(is.numeric(mfrow), length(mfrow) == 2)
nrow = mfrow[1]
ncol = mfrow[2]
}
xsize = sz[1, nrow]
ysize = sz[2, nrow]
asp = ysize / xsize
ret = list()
ret$mfrow = c(nrow, ncol)
# the following is right now only used by stars; FIXME:
# nocov start
ret$key.pos = if (!is.null(key.pos) && key.pos == -1L) { # figure out here: right or bottom?
newasp = asp * nrow / ncol # of the composition
dispasp = total_size[2] / total_size[1]
ifelse(newasp > dispasp, 4, 1)
} else
key.pos
m = matrix(seq_len(nrow * ncol), nrow, ncol, byrow = TRUE)
if (!is.null(ret$key.pos) && ret$key.pos != 0) { # add key row or column:
k = key.width
n = nrow * ncol + 1
switch(ret$key.pos,
{ ret$m = rbind(m, n); ret$widths = c(rep(1, ncol)); ret$heights = c(rep(asp, nrow), k) },
{ ret$m = cbind(n, m); ret$widths = c(k, rep(1, ncol)); ret$heights = c(rep(asp, nrow)) },
{ ret$m = rbind(n, m); ret$widths = c(rep(1, ncol)); ret$heights = c(k, rep(asp, nrow)) },
{ ret$m = cbind(m, n); ret$widths = c(rep(1, ncol), k); ret$heights = c(rep(asp, nrow)) }
)
} else {
ret$m = m
ret$widths = rep(1, ncol)
ret$heights = rep(asp, nrow)
}
# nocov end
ret
}
bb2merc = function(x, cls = "ggmap") { # return bbox in the appropriate "web mercator" CRS
wgs84 = st_crs(4326)
merc = st_crs(3857) # http://wiki.openstreetmap.org/wiki/EPSG:3857
pts = if (cls == "ggmap") {
b = vapply(attr(x, "bb"), c, 0.0)
st_sfc(st_point(c(b[2:1])), st_point(c(b[4:3])), crs = wgs84)
} else if (cls == "RgoogleMaps")
st_sfc(st_point(rev(x$BBOX$ll)), st_point(rev(x$BBOX$ur)), crs = wgs84)
else
stop("unknown cls")
st_bbox(st_transform(pts, merc))
}
#' @rdname stars
#' @param side ignore
#' @param at ignore
#' @param labels ignore
#' @param lon ignore
#' @param lat ignore
#' @param ndiscr ignore
#' @param reset ignore
#' @export
.degAxis = function (side, at, labels, ..., lon, lat, ndiscr, reset) {
if (missing(at))
at = axTicks(side)
if (missing(labels)) {
labels = FALSE
if (side == 1 || side == 3)
labels = parse(text = degreeLabelsEW(at))
else if (side == 2 || side == 4)
labels = parse(text = degreeLabelsNS(at))
}
axis(side, at = at, labels = labels, ...)
}
# find out where to place the legend key:
# given range r = (a, b), key.length l, key offset o, return a value range that:
# (i) scales such that (b - a) / (y - x) = l, and
# (ii) shifts linearly within [x, y] from a = x when o = 0 to b = y when o = 1
xy_from_r = function(r, l, o) {
stopifnot(length(r) == 2, l <= 1, l > 0, o >= 0, o <= 1)
r = as.numeric(r)
a = r[1]; b = r[2]
if (o == 1) {
y = b
x = b - (b - a)/l
} else {
i = o / (o - 1)
y = (a + (b - a)/l - i * b)/(1 - i)
x = i * (y - b) + a
}
c(x, y)
}
#' @rdname stars
#' @param z ignore
#' @param col ignore
#' @param breaks ignore
#' @param key.pos ignore
#' @param add.axis ignore
#' @param axes ignore
#' @param logz ignore
#' @param ... ignore
#' @param lab ignore
#' @param cex.axis see \link{par}
#' @export
.image_scale = function(z, col, breaks = NULL, key.pos, add.axis = TRUE,
at = NULL, ..., axes = FALSE, key.length, logz = FALSE, lab = "",
cex.axis = par("cex.axis")) {
if (!is.null(breaks) && length(breaks) != (length(col) + 1))
stop("must have one more break than colour")
stopifnot(is.null(lab) || is.character(lab) || is.expression(lab))
lab_set = (is.character(lab) && lab != "") || is.expression(lab)
zlim = range(z, na.rm = TRUE)
if (is.null(breaks))
breaks = seq(zlim[1], zlim[2], length.out = length(col) + 1)
offset = 0.5
if (length(key.pos) == 2) {
offset = key.pos[2]
key.pos = key.pos[1]
}
if (is.character(key.length)) {
kl = as.numeric(gsub(" cm", "", key.length))
sz = if (key.pos %in% c(1,3))
dev.size("cm")[1]
else
dev.size("cm")[2]
key.length = kl/sz
}
if (is.null(at)) {
br = range(breaks)
at = pretty(br)
at = at[at > br[1] & at < br[2]]
}
if (key.pos %in% c(1,3)) {
ylim = c(0, 1)
xlim = xy_from_r(range(breaks), key.length, offset)
mar = c(0, ifelse(axes, 2.1, 1), 0, 1)
}
if (key.pos %in% c(2,4)) {
ylim = xy_from_r(range(breaks), key.length, offset)
xlim = c(0, 1)
mar = c(ifelse(axes, 2.1, 1), 0, 1.2, 0)
}
mar[key.pos] = 2.1 + 1.5 * lab_set
par(mar = mar)
plot(1, 1, type = "n", ylim = ylim, xlim = xlim, axes = FALSE,
xlab = "", ylab = "", xaxs = "i", yaxs = "i")
if (!is.null(lab) && lab != "")
mtext(lab, side = key.pos, line = 2.5, cex = .8)
poly = vector(mode="list", length(col))
for (i in seq(poly))
poly[[i]] = c(breaks[i], breaks[i+1], breaks[i+1], breaks[i])
offset = 0.2
offs = switch(key.pos,
c(0, 0, -offset, -offset),
c(0, 0, -offset, -offset),
c(offset, offset, 0, 0),
c(offset, offset, 0, 0))
for(i in seq_along(poly)) {
if (key.pos %in% c(1,3))
polygon(poly[[i]], c(0, 0, 1, 1) + offs, col = col[i], border = NA)
if (key.pos %in% c(2,4))
polygon(c(0, 0, 1, 1) + offs, poly[[i]], col = col[i], border = NA)
}
# box() now would draw around [0,1]:
bx = c(breaks[1], rep(tail(breaks, 1), 2), breaks[1])
if (key.pos %in% c(1, 3))
polygon(bx, c(0, 0, 1, 1) + offs, col = NA, border = 'black')
if (key.pos %in% c(2, 4))
polygon(c(0, 0, 1, 1) + offs, bx, col = NA, border = 'black')
labels = if (logz)
parse(text = paste0("10^", at))
else if (inherits(breaks, c("POSIXt", "Date")))
format(at)
else
TRUE
if (add.axis)
axis(key.pos, at = at, labels = labels, cex.axis = cex.axis)
}
#' @rdname stars
#' @param key.width ignore
#' @export
.image_scale_factor = function(z, col, key.pos, add.axis = TRUE,
..., axes = FALSE, key.width, key.length, cex.axis = par("cex.axis")) {
to_num = function(x) as.numeric(gsub(" cm", "", x))
n = length(z)
if ((kw <- to_num(key.width)) < to_num(kw_dflt(list(z), key.pos))) { # cut z until it fits:
m = max(nchar(z))
while(kw < to_num(kw_dflt(list(z), key.pos))) { # cut z:
m = m - 1
z = substr(z, 1, m)
}
}
ksz = max(1.5 + max(nchar(z))/2, max(strwidth(z, "inches")) / par("cin")[1]) # in "mar" lines
breaks = (0:n) + 0.5
offset = 0.5
if (length(key.pos) == 2) {
offset = key.pos[2]
key.pos = key.pos[1]