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
}

Read our documentation on viewing source code .

Loading