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
|
R version 3.4.2 (2017-09-28) -- "Short Summer"
Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(maps)
> m = map('usa', plot = FALSE, fill = TRUE)
> suppressPackageStartupMessages(library(sf))
> m0 <- st_as_sfc(m)
> m <- st_as_sf(m)
>
> laea = st_crs("+proj=laea +lat_0=30 +lon_0=-95") # Lambert equal area
> m <- st_transform(st_as_sf(m), laea)
>
> bb = st_bbox(m)
> bbox = st_linestring(rbind(c( bb[1],bb[2]),c( bb[3],bb[2]),c( bb[3],bb[4]),c( bb[1],bb[4]),c( bb[1],bb[2])))
>
> g = st_graticule(m)
> plot(m, xlim = 1.2 * c(-2450853.4, 2186391.9))
> plot(g[1], add = TRUE, col = 'grey')
> plot(bbox, add = TRUE)
> points(g$x_start, g$y_start, col = 'red')
> points(g$x_end, g$y_end, col = 'blue')
>
> invisible(lapply(seq_len(nrow(g)), function(i) {
+ if (g$type[i] == "N" && g$x_start[i] - min(g$x_start) < 1000)
+ text(g[i,"x_start"], g[i,"y_start"], labels = parse(text = g[i,"degree_label"]),
+ srt = g$angle_start[i], pos = 2, cex = .7)
+ if (g$type[i] == "E" && g$y_start[i] - min(g$y_start) < 1000)
+ text(g[i,"x_start"], g[i,"y_start"], labels = parse(text = g[i,"degree_label"]),
+ srt = g$angle_start[i] - 90, pos = 1, cex = .7)
+ if (g$type[i] == "N" && g$x_end[i] - max(g$x_end) > -1000)
+ text(g[i,"x_end"], g[i,"y_end"], labels = parse(text = g[i,"degree_label"]),
+ srt = g$angle_end[i], pos = 4, cex = .7)
+ if (g$type[i] == "E" && g$y_end[i] - max(g$y_end) > -1000)
+ text(g[i,"x_end"], g[i,"y_end"], labels = parse(text = g[i,"degree_label"]),
+ srt = g$angle_end[i] - 90, pos = 3, cex = .7)
+ }))
>
> plot(m, graticule = st_crs(4326))
> nc = st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
> # options(warn=2)
> g = st_graticule(nc, datum = st_crs(nc))
> #g = st_graticule(nc)
>
> plot(nc[1], graticule = st_crs(nc))
>
> plot(nc[1], graticule = st_crs(nc), axes = TRUE)
>
> g = st_graticule()
>
> library(ggplot2)
> if (utils::packageVersion("ggplot2") > "2.2.1") {
+ ggplot() + geom_sf(data = st_set_crs(nc, NA_crs_)) # NA_crs_ for crs
+ }
>
> library(maps) #421
> wrld2 = st_as_sf(map('world2', plot=FALSE, fill=TRUE ))
> plot(wrld2, graticule = TRUE)
>
> proc.time()
user system elapsed
2.660 0.320 2.697
|
