|
1
|
|
rasterCheckSize <- function(x, interactive) {
|
|
2
|
|
# if (maxpixels < raster::ncell(x)) {
|
|
3
|
|
# warning(paste("maximum number of pixels for Raster* viewing is",
|
|
4
|
|
# maxpixels, "; \nthe supplied Raster* has", ncell(x), "\n",
|
|
5
|
|
# "... decreasing Raster* resolution to", maxpixels, "pixels\n",
|
|
6
|
|
# "to view full resolution set 'maxpixels = ", ncell(x), "'"))
|
|
7
|
|
|
|
8
|
2
|
tmapOptions <- get(".tmapOptions", envir = .TMAP_CACHE)
|
|
9
|
2
|
max.raster <- tmapOptions$max.raster
|
|
10
|
2
|
show.messages <- tmapOptions$show.messages
|
|
11
|
|
|
|
12
|
2
|
nc <- raster::ncell(x)
|
|
13
|
2
|
mx <- max.raster[ifelse(interactive, "view", "plot")]
|
|
14
|
2
|
if (nc > mx) {
|
|
15
|
0
|
if (show.messages) message("Raster object has ", nc, " (", nrow(x), " by ", ncol(x), ") cells, which is larger than ", mx, ", the maximum size determined by the option max.raster. Therefore, the raster will be shown at a decreased resolution of ", mx, " cells. Set tmap_options(max.raster = c(plot = ", nc, ", view = ", nc, ")) to show the whole raster.")
|
|
16
|
0
|
if (nlayers(x) > 1) {
|
|
17
|
0
|
x <- do.call(brick, lapply(1L:nlayers(x), function(i) {
|
|
18
|
0
|
raster::sampleRegular(raster(x, layer = i), mx, asRaster = TRUE, useGDAL = TRUE)
|
|
19
|
0
|
}))
|
|
20
|
0
|
} else {
|
|
21
|
0
|
x <- raster::sampleRegular(x, mx, asRaster = TRUE, useGDAL = TRUE)
|
|
22
|
0
|
}
|
|
23
|
2
|
}
|
|
24
|
|
|
|
25
|
|
|
|
26
|
2
|
return(x)
|
|
27
|
|
}
|
|
28
|
|
|
|
29
|
|
preprocess_shapes <- function(y, raster_facets_vars, gm, interactive) {
|
|
30
|
2
|
shp <- y$shp
|
|
31
|
|
|
|
32
|
2
|
show.messages <- get(".tmapOptions", envir = .TMAP_CACHE)$show.messages
|
|
33
|
|
|
|
34
|
|
|
|
35
|
2
|
if (is.null(shp)) return(list(shp=NULL, data=NULL, type="tiles"))
|
|
36
|
|
|
|
37
|
2
|
shp.unit <- gm$shape.unit
|
|
38
|
|
# shp.aa <- y[c("unit", "orig", "to", "total.area")]
|
|
39
|
|
# names(shp.aa)[names(shp.aa)=="unit"] <- "target"
|
|
40
|
|
# shp.aa <- shp.aa[!sapply(shp.aa, is.null)]
|
|
41
|
|
|
|
42
|
2
|
shp.sim <- y[c("simplify", "keep.units", "keep.subunits", "method", "no_repair", "snap", "force_FC", "drop_null_geometries")]
|
|
43
|
2
|
names(shp.sim)[names(shp.sim)=="simplify"] <- "fact"
|
|
44
|
2
|
shp.sim <- shp.sim[!vapply(shp.sim, is.null, logical(1))]
|
|
45
|
|
|
|
46
|
2
|
if (inherits(shp, c("Raster", "SpatialPixels", "SpatialGrid"))) {
|
|
47
|
2
|
is.RGB <- attr(raster_facets_vars, "is.RGB") # true if tm_rgb is used (NA if qtm is used)
|
|
48
|
2
|
rgb.vars <- attr(raster_facets_vars, "rgb.vars")
|
|
49
|
2
|
to.Cat <- attr(raster_facets_vars, "to.Cat") # true if tm_raster(..., style = "cat) is specified
|
|
50
|
2
|
max.value <- attr(raster_facets_vars, "max.value") # NULL is tm_raster is called, when tm_rgb is called: NA (default) when max color value is determined automatically.
|
|
51
|
|
|
|
52
|
2
|
if (interactive) gm$shape.master_crs <- .crs_merc
|
|
53
|
|
|
|
54
|
2
|
if (inherits(shp, "Spatial")) shp <- brick(shp)
|
|
55
|
|
|
|
56
|
|
|
|
57
|
|
|
|
58
|
|
# attribute get from read_osm
|
|
59
|
2
|
is.OSM <- attr(shp, "is.OSM")
|
|
60
|
2
|
if (is.null(is.OSM)) is.OSM <- FALSE
|
|
61
|
2
|
leaflet.server <- attr(shp, "leaflet.provider")
|
|
62
|
2
|
if (is.null(leaflet.server)) leaflet.server <- NA
|
|
63
|
|
|
|
64
|
|
# color values are encoded by a colortable (and not interpreted as factors)
|
|
65
|
2
|
if (length(colortable(shp))>0) {
|
|
66
|
0
|
ctable <- colortable(shp)
|
|
67
|
0
|
uctable <- unique(ctable)
|
|
68
|
0
|
mtch <- match(ctable, uctable)
|
|
69
|
|
|
|
70
|
0
|
if (nlayers(shp)>1) shp <- raster::subset(shp, 1)
|
|
71
|
0
|
shp <- setValues(shp, mtch[getValues(shp) + 1L])
|
|
72
|
0
|
names(shp) <- shpnames <- "PIXEL__COLOR"
|
|
73
|
|
|
|
74
|
|
|
|
75
|
0
|
use_interp <- FALSE
|
|
76
|
|
|
|
77
|
0
|
lvls <- list(uctable)
|
|
78
|
|
|
|
79
|
|
# if (!is.RGB && is.na(do.interpolate)) {
|
|
80
|
|
# if (get(".tmapOptions", envir = .TMAP_CACHE)$show.messages) {
|
|
81
|
|
# message("For bitmap images, it is recommended to use tm_rgb instead of tm_raster (or to set interpolate to TRUE).")
|
|
82
|
|
# }
|
|
83
|
|
# }
|
|
84
|
0
|
layerIDs <- 1
|
|
85
|
0
|
convert.RGB <- FALSE
|
|
86
|
2
|
} else {
|
|
87
|
|
# in order to not loose factor levels, subset the data here
|
|
88
|
2
|
rdata <- get_raster_data(shp, show.warnings = FALSE)
|
|
89
|
2
|
mainID <- attr(rdata, "mainID")
|
|
90
|
|
|
|
91
|
|
|
|
92
|
2
|
shpnames <- names(rdata) #get_raster_names(shp)
|
|
93
|
|
|
|
94
|
2
|
mxdata <- max(maxValue(shp))
|
|
95
|
|
|
|
96
|
2
|
if (is.na(is.RGB)) {
|
|
97
|
2
|
if ((nlayers(shp) == 3) && all(minValue(shp)>=0) && mxdata <= max.value) {
|
|
98
|
0
|
if (mxdata <= 1) {
|
|
99
|
0
|
max.value <- 1
|
|
100
|
0
|
message("Numeric values of ", y$shp_name, " interpreted as RGB values with max.value = 1. Run tm_shape(", y$shp_name, ") + tm_raster() to visualize the data.")
|
|
101
|
0
|
} else {
|
|
102
|
0
|
message("Numeric values of ", y$shp_name, " interpreted as RGB values with max.value = 255. Run tm_shape(", y$shp_name, ") + tm_raster() to visualize the data.")
|
|
103
|
0
|
}
|
|
104
|
0
|
is.RGB <- TRUE
|
|
105
|
0
|
rgb.vars <- 1:3
|
|
106
|
2
|
} else {
|
|
107
|
2
|
is.RGB <- FALSE
|
|
108
|
2
|
}
|
|
109
|
2
|
} else if (is.RGB) {
|
|
110
|
0
|
if (!any(rgb.vars %in% 1:nlayers(shp))) stop("Specified rgb(a) bands are ", rgb.vars, " whereas the number of layers is ", nlayers(shp), call. = FALSE)
|
|
111
|
0
|
if (!all(minValue(shp)>=0) || !all(mxdata <= max.value)) {
|
|
112
|
0
|
shp[][shp[] < 0] <- 0
|
|
113
|
0
|
shp[][shp[] > max.value] <- max.value
|
|
114
|
0
|
rdata[rdata < 0] <- 0
|
|
115
|
0
|
rdata[rdata > max.value] <- max.value
|
|
116
|
0
|
warning("Raster values found that are outside the range [0, ", max.value, "]", call. = FALSE)
|
|
117
|
0
|
}
|
|
118
|
0
|
if (mxdata <= 1 && max.value == 255) message("No values higher than 1 found. Probably, max.value should be set to 1.")
|
|
119
|
2
|
}
|
|
120
|
|
|
|
121
|
2
|
convert.RGB <- is.RGB
|
|
122
|
|
|
|
123
|
|
|
|
124
|
|
#
|
|
125
|
|
#
|
|
126
|
|
# convert.RGB <- !identical(is.RGB, FALSE) &&
|
|
127
|
|
# (is.na(raster_facets_vars[1]) || !any(raster_facets_vars %in% shpnames)) &&
|
|
128
|
|
# all(minValue(shp)>=0) && mxdata <= max.value
|
|
129
|
|
#
|
|
130
|
|
# if (is.na(is.RGB) && convert.RGB && show.messages) {
|
|
131
|
|
# message("Numeric values of ", y$shp_name, " interpreted as RGB(A) values. Run tm_shape(", y$shp_name, ") + tm_raster() to visualize the data.")
|
|
132
|
|
# }
|
|
133
|
|
#
|
|
134
|
|
#
|
|
135
|
|
# if (identical(is.RGB, TRUE) && !convert.RGB) {
|
|
136
|
|
# stop("Raster object does not have a color table, nor numeric data that can be converted to colors. Use tm_raster to visualize the data.", call. = FALSE)
|
|
137
|
|
# }
|
|
138
|
|
|
|
139
|
2
|
if (convert.RGB) {
|
|
140
|
0
|
layerIDs <- rgb.vars #1L:nlayers(shp)
|
|
141
|
2
|
} else {
|
|
142
|
2
|
if (is.na(raster_facets_vars[1])) {
|
|
143
|
2
|
if (length(mainID) != length(shpnames) && show.messages) {
|
|
144
|
0
|
if (attr(rdata, "cls") == "RasterLayer") {
|
|
145
|
0
|
message("Only the first variable is shown. The available variables are: \"", paste(shpnames, collapse = "\", \""), "\".")
|
|
146
|
0
|
} else {
|
|
147
|
0
|
message("For each raster layer, only the first variable is shown. The available variables are: \"", paste(shpnames, collapse = "\", \""), "\".")
|
|
148
|
0
|
}
|
|
149
|
2
|
}
|
|
150
|
2
|
raster_facets_vars <- shpnames[mainID]
|
|
151
|
2
|
} else {
|
|
152
|
0
|
raster_facets_vars <- na.omit(raster_facets_vars)
|
|
153
|
2
|
}
|
|
154
|
|
|
|
155
|
2
|
layerIDs <- match(raster_facets_vars, shpnames)
|
|
156
|
2
|
layerIDs <- na.omit(layerIDs)
|
|
157
|
2
|
if (length(layerIDs) == 0L) layerIDs <- 1
|
|
158
|
2
|
}
|
|
159
|
|
#lvls <- get_raster_levels(shp, layerIDs)
|
|
160
|
2
|
lvls <- get_data_frame_levels(rdata[, layerIDs, drop = FALSE])
|
|
161
|
|
|
|
162
|
2
|
use_interp <- ((all(vapply(lvls, is.null, logical(1)))) && !to.Cat)
|
|
163
|
2
|
}
|
|
164
|
|
|
|
165
|
|
#print(shp)
|
|
166
|
|
|
|
167
|
2
|
shp <- rasterCheckSize(shp, interactive)
|
|
168
|
|
|
|
169
|
|
# print(use_interp)
|
|
170
|
|
# print(lvls)
|
|
171
|
|
# print(layerIDs)
|
|
172
|
|
# print(use_interp)
|
|
173
|
|
#shp <- rasterCheckSize(shp)
|
|
174
|
|
|
|
175
|
|
|
|
176
|
|
# get current projection (assume longlat if unkown)
|
|
177
|
2
|
shp_crs <- get_projection(shp, output="crs")
|
|
178
|
2
|
if (is.na(shp_crs)) {
|
|
179
|
0
|
if (!tmaptools::is_projected(shp)) {
|
|
180
|
0
|
warning("Currect projection of shape ", y$shp_name, " unknown. Long-lat (WGS84) is assumed.", call. = FALSE)
|
|
181
|
0
|
shp_crs <- .crs_longlat
|
|
182
|
0
|
shp <- set_projection(shp, current.projection = shp_crs)
|
|
183
|
0
|
} else {
|
|
184
|
0
|
warning("Current projection of shape ", y$shp_name, " unknown and cannot be determined.", call. = FALSE)
|
|
185
|
0
|
}
|
|
186
|
2
|
}
|
|
187
|
|
|
|
188
|
|
# should raster shape be reprojected?
|
|
189
|
2
|
if ((!is.na(shp_crs) && !is.na(gm$shape.master_crs) && !identical(shp_crs$proj4string, gm$shape.master_crs$proj4string)) || interactive) {
|
|
190
|
0
|
if (is.na(gm$shape.master_crs)) stop("Master projection unknown, but needed to reproject raster shape.", call.=FALSE)
|
|
191
|
0
|
new_ext <- tryCatch({
|
|
192
|
0
|
suppressWarnings(projectExtent(shp, crs = gm$shape.master_crs$proj4string))
|
|
193
|
0
|
}, error=function(e){
|
|
194
|
0
|
NULL
|
|
195
|
0
|
})
|
|
196
|
|
|
|
197
|
0
|
if (is.null(new_ext)) {
|
|
198
|
0
|
shp <- crop_shape(shp, bb(gm$shape.bbx_raw, projection = shp_crs, current.projection = gm$shape.master_crs))
|
|
199
|
0
|
new_ext <- tryCatch({
|
|
200
|
0
|
suppressWarnings(projectExtent(shp, crs = gm$shape.master_crs$proj4string))
|
|
201
|
0
|
}, error=function(e){
|
|
202
|
0
|
stop("Unable to reproject raster shape \"", y$shp_name, "\", probably due to non-finite points.", call. = FALSE)
|
|
203
|
0
|
})
|
|
204
|
0
|
}
|
|
205
|
2
|
} else new_ext <- NULL
|
|
206
|
|
|
|
207
|
|
|
|
208
|
2
|
raster.projected <- !is.null(new_ext)
|
|
209
|
|
|
|
210
|
2
|
if (raster.projected) {
|
|
211
|
0
|
shpTmp <- suppressWarnings(projectRaster(shp, to=new_ext, crs=gm$shape.master_crs$proj4string, method = ifelse(use_interp, "bilinear", "ngb")))
|
|
212
|
0
|
shp2 <- raster(shpTmp)
|
|
213
|
0
|
data <- suppressWarnings(get_raster_data(shpTmp)[,layerIDs, drop=FALSE])
|
|
214
|
2
|
} else {
|
|
215
|
2
|
shp2 <- raster(shp)
|
|
216
|
2
|
data <- suppressWarnings(get_raster_data(shp)[,layerIDs, drop=FALSE])
|
|
217
|
2
|
}
|
|
218
|
|
|
|
219
|
|
# restore factor levels and limits
|
|
220
|
2
|
data <- as.data.frame(mapply(function(d, l) {
|
|
221
|
2
|
if (is.character(l) && !is.factor(d)) {
|
|
222
|
0
|
if (is.logical(d)) {
|
|
223
|
0
|
d2 <- as.integer(d)+1L
|
|
224
|
0
|
} else {
|
|
225
|
0
|
plusone <- min(d, na.rm=TRUE)==0
|
|
226
|
0
|
if (!is.integer(d)) {
|
|
227
|
0
|
d2 <- as.integer(d) + plusone
|
|
228
|
0
|
} else {
|
|
229
|
0
|
d2 <- d+plusone
|
|
230
|
0
|
}
|
|
231
|
0
|
}
|
|
232
|
0
|
levels(d2) <- l
|
|
233
|
0
|
class(d2) <- "factor"
|
|
234
|
0
|
d2
|
|
235
|
2
|
} else if (is.numeric(l)) {
|
|
236
|
2
|
pmin(pmax(l[1], d), l[2])
|
|
237
|
2
|
} else d
|
|
238
|
2
|
}, data, lvls, SIMPLIFY=FALSE))
|
|
239
|
|
|
|
240
|
2
|
if (convert.RGB) {
|
|
241
|
0
|
data <- data.frame(PIXEL__COLOR = raster_colors(as.matrix(data), use.colortable = FALSE, max.value = max.value))
|
|
242
|
2
|
}
|
|
243
|
|
|
|
244
|
|
|
|
245
|
|
# set values in order to align data later on (with cropping)
|
|
246
|
2
|
shp2 <- setValues(shp2, values=1:ncell(shp2))
|
|
247
|
|
|
|
248
|
|
# interactive raster require merc projection with longlat extent
|
|
249
|
2
|
if (interactive) {
|
|
250
|
0
|
new_ext_ll <- extent(projectExtent(new_ext, crs = .crs_longlat$proj4string))
|
|
251
|
0
|
shp2@extent <- new_ext_ll
|
|
252
|
0
|
shp2@crs <- raster::crs(.crs_longlat$proj4string)
|
|
253
|
2
|
}
|
|
254
|
|
|
|
255
|
|
## to be consistent with Spatial objects:
|
|
256
|
2
|
attr(shp2, "bbox") <- bb(shp2)
|
|
257
|
2
|
attr(shp2, "proj4string") <- attr(shp2@crs, "projargs")
|
|
258
|
|
|
|
259
|
|
|
|
260
|
2
|
data$tmapfilter <- TRUE
|
|
261
|
|
|
|
262
|
2
|
attr(data, "is.OSM") <- is.OSM
|
|
263
|
2
|
attr(data, "leaflet.server") <- leaflet.server
|
|
264
|
2
|
attr(data, "raster.projected") <- raster.projected
|
|
265
|
|
|
|
266
|
|
#attr(data, "is.RGB") <- is.RGB
|
|
267
|
|
|
|
268
|
2
|
type <- "raster"
|
|
269
|
|
|
|
270
|
2
|
} else {
|
|
271
|
|
# save_bbox (sp objects allow for custom bboxes, sf objects don't)
|
|
272
|
2
|
shp_bbx <- bb(shp)
|
|
273
|
|
|
|
274
|
2
|
kernel_density <- ("kernel_density" %in% names(attributes(shp)))
|
|
275
|
2
|
isolines <- ("isolines" %in% names(attributes(shp)))
|
|
276
|
|
|
|
277
|
2
|
if (inherits(shp, "Spatial")) {
|
|
278
|
0
|
shp <- as(shp, "sf")
|
|
279
|
2
|
} else if (!inherits(shp, c("sf", "sfc"))) {
|
|
280
|
0
|
stop("Object ", y$shp_name, " is neither from class sf, Spatial, nor Raster.", call. = FALSE)
|
|
281
|
2
|
}
|
|
282
|
|
|
|
283
|
|
# drop z/m
|
|
284
|
2
|
shp <- sf::st_zm(shp)
|
|
285
|
|
|
|
286
|
|
# remove empty units
|
|
287
|
2
|
empty_units <- st_is_empty(shp)
|
|
288
|
2
|
if (any(empty_units)) {
|
|
289
|
0
|
shp <- if (inherits(shp, "sf")) shp[!empty_units, ] else shp[!empty_units]
|
|
290
|
2
|
}
|
|
291
|
|
|
|
292
|
|
|
|
293
|
|
## get data.frame from shapes, and store ID numbers in shape objects (needed for cropping)
|
|
294
|
2
|
if (inherits(shp, "sfc")) {
|
|
295
|
0
|
data <- data.frame(tmapID = seq_len(length(shp)))
|
|
296
|
0
|
if (!is.null(names(shp))) names(shp) <- NULL
|
|
297
|
0
|
shp <- st_sf(data, geometry=shp)
|
|
298
|
2
|
} else {
|
|
299
|
2
|
data <- shp
|
|
300
|
2
|
st_geometry(data) <- NULL
|
|
301
|
2
|
shp <- st_geometry(shp)
|
|
302
|
2
|
if (!is.null(names(shp))) names(shp) <- NULL
|
|
303
|
2
|
shp <- st_sf(tmapID = seq_len(length(shp)), geometry = shp)
|
|
304
|
2
|
}
|
|
305
|
|
|
|
306
|
2
|
data$tmapfilter <- if (is.null(y$filter)) rep(TRUE, nrow(shp)) else rep(y$filter, length.out = nrow(shp))
|
|
307
|
|
|
|
308
|
|
# reproject if nessesary
|
|
309
|
2
|
shp_crs <- get_projection(shp, output="crs")
|
|
310
|
2
|
if (is.na(shp_crs)) {
|
|
311
|
0
|
if (!tmaptools::is_projected(shp)) {
|
|
312
|
0
|
warning("Currect projection of shape ", y$shp_name, " unknown. Long-lat (WGS84) is assumed.", call. = FALSE)
|
|
313
|
0
|
shp_crs <- .crs_longlat
|
|
314
|
0
|
shp <- set_projection(shp, current.projection = shp_crs)
|
|
315
|
0
|
} else {
|
|
316
|
0
|
warning("Current projection of shape ", y$shp_name, " unknown and cannot be determined.", call. = FALSE)
|
|
317
|
0
|
}
|
|
318
|
2
|
}
|
|
319
|
2
|
if (!is.na(shp_crs) && !is.na(gm$shape.master_crs) && !identical(shp_crs$proj4string, gm$shape.master_crs$proj4string)) {
|
|
320
|
2
|
shp2 <- set_projection(shp, gm$shape.master_crs)
|
|
321
|
|
|
|
322
|
|
# override bounding box (since it now is projected)
|
|
323
|
2
|
shp_bbx <- bb(shp2)
|
|
324
|
2
|
} else {
|
|
325
|
2
|
shp2 <- shp
|
|
326
|
2
|
}
|
|
327
|
|
|
|
328
|
2
|
if (inherits(st_geometry(shp2), c("sfc_POLYGON", "sfc_MULTIPOLYGON"))) {
|
|
329
|
2
|
data$SHAPE_AREAS <- tmaptools::approx_areas(shp=shp2, target = paste(shp.unit, shp.unit, sep=" "))
|
|
330
|
2
|
if (gm$shape.apply_map_coloring) attr(data, "NB") <- if (length(shp)==1) list(0) else get_neighbours(shp) #poly2nb(as(shp, "Spatial"))
|
|
331
|
2
|
attr(data, "kernel_density") <- kernel_density
|
|
332
|
2
|
type <- "polygons"
|
|
333
|
2
|
} else if (inherits(st_geometry(shp2), c("sfc_LINESTRING", "sfc_MULTILINESTRING"))) {
|
|
334
|
2
|
attr(data, "isolines") <- isolines
|
|
335
|
|
## TODO update smooth_map to sf
|
|
336
|
2
|
type <- "lines"
|
|
337
|
2
|
} else if (inherits(st_geometry(shp2), c("sfc_POINT", "sfc_MULTIPOINT"))){
|
|
338
|
2
|
type <- "points"
|
|
339
|
2
|
} else {
|
|
340
|
2
|
if (any(st_geometry_type(shp2) == "GEOMETRYCOLLECTION")) {
|
|
341
|
0
|
gnew <- split_geometry_collection(st_geometry(shp2))
|
|
342
|
0
|
ids <- attr(gnew, "ids")
|
|
343
|
0
|
data <- data[ids, , drop = FALSE]
|
|
344
|
0
|
shp2 <- st_sf(tmapID = 1L:nrow(data), geometry = gnew)
|
|
345
|
2
|
}
|
|
346
|
2
|
type <- "geometrycollection"
|
|
347
|
2
|
attr(data, "kernel_density") <- FALSE
|
|
348
|
2
|
attr(type, "types") <- get_types(st_geometry(shp2))
|
|
349
|
2
|
}
|
|
350
|
|
|
|
351
|
|
# simplify shape
|
|
352
|
|
|
|
353
|
2
|
if (shp.sim$fact != 1 && type %in% c("polygons", "lines")) {
|
|
354
|
|
## TODO convert fact to tolerance
|
|
355
|
|
|
|
356
|
0
|
if (!requireNamespace("rmapshaper", quietly = TRUE)) {
|
|
357
|
0
|
warning("rmapshaper package is needed to simplify the shape. Alternatively, st_simplify from the sf package can be used. See the underlying function tmaptools::simplify_shape for details.", call. = FALSE)
|
|
358
|
0
|
} else {
|
|
359
|
|
#shp2 <- st_simplify(shp2, preserveTopology = TRUE, dTolerance = shp.sim$fact)
|
|
360
|
0
|
shp2 <- do.call(tmaptools::simplify_shape, c(list(shp=shp2), shp.sim))
|
|
361
|
0
|
data <- data[shp2$tmapID, , drop=FALSE]
|
|
362
|
0
|
shp2$tmapID <- seq_len(nrow(shp2))
|
|
363
|
0
|
}
|
|
364
|
2
|
}
|
|
365
|
|
|
|
366
|
|
# be consistent with rasters (originated from sp objects)
|
|
367
|
2
|
attr(shp2, "bbox") <- shp_bbx
|
|
368
|
2
|
attr(shp2, "proj4string") <- st_crs(shp2)
|
|
369
|
|
|
|
370
|
2
|
shpnames <- names(data)
|
|
371
|
|
|
|
372
|
2
|
}
|
|
373
|
|
|
|
374
|
2
|
point.per <- if (is.na(y$point.per)) ifelse(type %in% c("points", "geometrycollection"), "segment", "feature") else y$point.per
|
|
375
|
|
|
|
376
|
2
|
attr(data, "shpnames") <- shpnames
|
|
377
|
|
|
|
378
|
2
|
attr(shp2, "point.per") <- point.per
|
|
379
|
2
|
attr(shp2, "line.center") <- y$line.center
|
|
380
|
2
|
attr(shp2, "projected") <- tmaptools::is_projected(shp2)
|
|
381
|
2
|
list(shp=shp2, data=data, type=type)
|
|
382
|
|
}
|
|
383
|
|
|
|
384
|
|
get_types <- function(sfc) {
|
|
385
|
2
|
tp <- st_geometry_type(sfc)
|
|
386
|
2
|
types <- factor(rep(NA, length(sfc)), levels=c("polygons", "lines", "points", "collection"))
|
|
387
|
2
|
types[tp %in% c("MULTIPOLYGON", "POLYGON")] <- "polygons"
|
|
388
|
2
|
types[tp %in% c("MULTILINESTRING", "LINESTRING")] <- "lines"
|
|
389
|
2
|
types[tp %in% c("MULTIPOINT", "POINT")] <- "points"
|
|
390
|
2
|
types[tp == "GEOMETRYCOLLECTION"] <- "collection"
|
|
391
|
2
|
if (any(is.na(types))) stop("The following geometry types are not supported: ", paste(unique(tp[is.na(types)]), collapse = ", "), call. = FALSE)
|
|
392
|
2
|
types
|
|
393
|
|
}
|
|
394
|
|
|
|
395
|
|
split_geometry_collection <- function(sfc) {
|
|
396
|
0
|
types <- get_types(sfc)
|
|
397
|
0
|
res <- mapply(function(g, tp, id) {
|
|
398
|
0
|
if (tp == "collection") {
|
|
399
|
0
|
g2 <- suppressWarnings(list(st_collection_extract(g, "POLYGON"),
|
|
400
|
0
|
st_collection_extract(g, "POINT"),
|
|
401
|
0
|
st_collection_extract(g, "LINESTRING")))
|
|
402
|
|
# tp2 <- factor(c("polygons", "points", "lines"), levels=c("polygons", "lines", "points"))
|
|
403
|
0
|
sel <- !vapply(g2, st_is_empty, logical(1))
|
|
404
|
|
|
|
405
|
0
|
g2 <- g2[sel]
|
|
406
|
|
# tp2 <- tp2[sel]
|
|
407
|
0
|
id2 <- rep(id, length(g2))
|
|
408
|
0
|
list(g2, id2)
|
|
409
|
0
|
} else {
|
|
410
|
0
|
list(list(g), id)
|
|
411
|
0
|
}
|
|
412
|
0
|
}, sfc, types, 1:length(sfc), SIMPLIFY = FALSE)
|
|
413
|
0
|
gnew <- st_sfc(do.call(c, lapply(res, "[[", 1)), crs = st_crs(sfc))
|
|
414
|
0
|
ids <- do.call(c, lapply(res, "[[", 2))
|
|
415
|
0
|
attr(gnew, "ids") <- ids
|
|
416
|
0
|
gnew
|
|
417
|
|
}
|
mtennekes