reorganized scripts, now split by spatial package

R/auto_crs.R

@@ -2,5 +2,6 @@ auto_crs = function(x, world_crs = "+proj=eck4") {
 	# placeholder for new automatic-crs function implemented elsewhere
 
 	# for the time being: use unprojected coordinates
-	4326
+	st_crs(x) # was 4326, but this may cause transformation
+
 }

R/spatial_raster.R

@@ -0,0 +1,14 @@
+#' @export
+tmapShape.Raster = function(shp, is.main, crs, bbox, unit, filter, shp_name, smeta, o, tmf) {
+	tmapShape.SpatRaster(terra::rast(shp), is.main, crs, bbox, unit, filter, shp_name, smeta, o, tmf)
+}
+
+#' @export
+tmapSubsetShp.Raster = function(shp, vars) {
+	tmapSubsetShp.SpatRaster(terra::rast(shp), vars)
+}
+
+#' @export
+tmapGetShapeMeta1.Raster = function(shp, o) {
+	tmapGetShapeMeta1.SpatRaster(terra::rast(shp), o)
+}

R/spatial_sf.R

@@ -0,0 +1,144 @@
+#' @export
+tmapReproject.sfc = function(shp, tmapID, bbox = NULL, ..., crs) {
+	if (is.na(sf::st_crs(shp))) {
+		shp2 = shp
+		#sf::st_crs(shp2) = crs
+		#warning("Setting missing CRS to ", as.character(crs))
+	} else {
+		shp2 = sf::st_transform(shp, crs)
+	}
+	if (!is.null(bbox$x)) bbox = list(x = do.call(tmaptools::bb, c(bbox, list(projection = crs))))
+	shapeTM(shp2, tmapID, bbox, ...)
+}
+
+#' @export
+tmapShape.sf = function(shp, is.main, crs, bbox, unit, filter, shp_name, smeta, o, tmf) {
+	if (identical(crs, "auto")) crs = auto_crs(shp) else crs = sf::st_crs(crs)
+	reproj = (!is.null(crs) && !is.na(crs) && sf::st_crs(shp) != crs)
+
+	if (reproj) {
+		if (crs_is_ortho(crs)) {
+			tryCatch({
+				suppressWarnings({
+					shp4326 = sf::st_transform(shp, 4326)
+					visible = crs_ortho_visible(crs, projected = FALSE)
+					if (!sf::st_is_valid(visible)) visible = sf::st_make_valid(visible)
+					shp = suppressMessages(sf::st_intersection(shp4326, visible))
+				})
+			}, error = function(e) {
+				shp
+			})
+		}
+		shp = sf::st_transform(shp, crs = crs)
+	}
+
+	sfc = sf::st_geometry(shp)
+
+	if (inherits(sfc[[1]], c("XYZ", "XYM", "XYZM"))) {
+		sfc = sf::st_zm(sfc)
+	}
+
+	if (o$check.and.fix) sfc = check_fix(sfc, shp_name, reproj, o$show.messages)
+
+	# if check_fix fails, is_valid contains the valid ids
+	isv = attr(sfc, "is_valid")
+	if (!is.null(isv)) {
+		shp = shp[isv,]
+	}
+	dt = data.table::as.data.table(sf::st_drop_geometry(shp))
+
+	dtcols = copy(names(dt))
+
+	#if (is.null(bbox)) bbox = sf::st_bbox(sfc)
+
+	if (is.null(filter)) filter = rep_len(TRUE, nrow(dt))
+	dt[, ':='(tmapID__ = 1L:nrow(dt), sel__ = filter)]
+
+	make_by_vars(dt, tmf, smeta)
+
+	shpTM = shapeTM(shp = sfc, tmapID = 1L:(length(sfc)), bbox = bbox)
+
+	structure(list(shpTM = shpTM, dt = dt, is.main = is.main, dtcols = dtcols, shpclass = "sfc", bbox = bbox, unit = unit, shp_name = shp_name, smeta = smeta), class = "tmapShape")
+}
+
+#' @export
+tmapSubsetShp.sf = function(shp, vars) {
+
+	if ("AREA" %in% vars && !("AREA" %in% names(shp))) {
+		shp$AREA = sf::st_area(shp)
+	}
+	if ("LENGTH" %in% vars && !("LENGTH" %in% names(shp))) {
+		shp$LENGTH = sf::st_length(shp)
+	}
+	if ("MAP_COLORS" %in% vars) {
+		shp$MAP_COLORS = as.factor(tmaptools::map_coloring(shp))
+	}
+
+	if (!length(vars)) {
+		vars = "dummy__"
+		shp$dummy__ = logical(nrow(shp))
+	}
+	shp[, vars]
+}
+
+
+#' @export
+tmapSubsetShp.sfc = function(shp, vars) {
+	s = sf::st_sf(dummy__ = TRUE, geometry = shp)
+	if ("AREA" %in% vars) {
+		s$AREA = sf::st_area(shp)
+	}
+	if ("LENGTH" %in% vars) {
+		s$LENGTH = sf::st_length(shp)
+	}
+	if ("MAP_COLORS" %in% vars) {
+		s$MAP_COLORS = tmaptools::map_coloring(s)
+	}
+	s
+}
+
+
+
+#' @export
+tmapGetShapeMeta2.sf = function(shp, smeta, o) {
+	vars = setdiff(names(shp), attr(shp, "sf_column"))
+	smeta$vars_levs = lapply(seq_along(vars), function(i) {
+		get_fact_levels_na(shp[[i]], o)
+	})
+	names(smeta$vars_levs) = vars
+	smeta
+}
+
+
+#' @export
+tmapGetShapeMeta1.sf = function(shp, o) {
+	vars = setdiff(names(shp), attr(shp, "sf_column"))
+	names(vars) = vars
+	#vars_levs = lapply(vars, function(v) {get_fact_levels_na(shp[[v]], o)})
+	dims = character(0)
+	dims_vals = list()
+
+	list(vars = vars,
+		 dims = dims,
+		 dims_vals = dims_vals)
+}
+
+
+#' @export
+tmapGetShapeMeta2.sfc = function(shp, smeta, o) {
+	vars = character(0)
+	smeta$vars_levs = list()
+	smeta
+}
+
+
+#' @export
+tmapGetShapeMeta1.sfc = function(shp, o) {
+	vars = character(0)
+	dims = character(0)
+	dims_vals = list()
+
+	list(vars = vars,
+		 dims = dims,
+		 dims_vals = dims_vals)
+}

R/spatial_sp.R

@@ -0,0 +1,15 @@
+#' @export
+tmapShape.Spatial = function(shp, is.main, crs, bbox, unit, filter, shp_name, smeta, o, tmf) {
+	tmapShape.sf(as(shp, "sf"), is.main, crs, bbox, unit, filter, shp_name, smeta, o, tmf)
+}
+
+#' @export
+tmapSubsetShp.Spatial = function(shp, vars) {
+	tmapSubsetShp.sf(as(shp, "sf"), vars)
+}
+
+
+#' @export
+tmapGetShapeMeta1.Spatial = function(shp, o) {
+	tmapGetShapeMeta1.SpatRaster(as(shp, "sf"), o)
+}

R/spatial_stars.R

@@ -0,0 +1,200 @@
+#' @export
+tmapReproject.stars = function(shp, tmapID, bbox = NULL, ..., crs) {
+	shp[[1]] = tmapID
+
+	shp2 = transwarp(shp, crs, raster.warp = TRUE)
+	tmapID2 = shp2[[1]]
+	shp2[[1]] = NA
+
+	if (!is.null(bbox$x)) bbox = list(x = do.call(tmaptools::bb, c(bbox, list(projection = crs))))
+
+	shapeTM(shp2, tmapID2, bbox, ...)
+}
+
+#' @export
+tmapReproject.dimensions = function(shp, tmapID, bbox = NULL, ..., crs) {
+	s = structure(list(tmapID = matrix(tmapID, nrow = nrow(shp))), dimensions = shp, class = "stars")
+
+	if (is.na(sf::st_crs(shp))) {
+		shp2 = s
+	} else {
+		shp2 = transwarp(s, crs, raster.warp = TRUE)
+	}
+
+	tmapID2 = shp2[[1]]
+
+	d2 = stars::st_dimensions(shp2)
+	if (!is.null(bbox$x)) bbox = list(x = do.call(tmaptools::bb, c(bbox, list(projection = crs))))
+
+	shapeTM(d2, tmapID2, bbox, ...)
+}
+
+
+#' @export
+tmapShape.stars = function(shp, is.main, crs, bbox, unit, filter, shp_name, smeta, o, tmf) {
+	if (identical(crs, "auto")) crs = auto_crs(shp)
+
+	dev = getOption("tmap.devel.mode")
+
+	if (dev) cat("-- stars object:", shp_name, "--\n")
+
+
+	if (!has_raster(shp)) {
+		dimnms = dimnames(shp)
+
+		dimvals = lapply(1:length(dimnms), function(i) stars::st_get_dimension_values(shp, i))
+		dimsfc = vapply(dimvals, inherits, what = "sfc", FUN.VALUE = logical(1))
+
+		if (!any(dimsfc)) {
+			stop("stars object ", shp_name, " is a stars object without raster and doens't have a geometry dimension")
+		} else {
+			dimid = which(dimsfc)
+			geoms = dimvals[[dimid]]
+			dimnms_new = dimnms
+			dimnms_new[dimid] = "tmapID__"
+			dimcols = dimnms_new[-dimid] # columns names, used for default facetting
+			dimvls = lapply(dimcols, function(d) stars::st_get_dimension_values(shp, d))
+			shpnames = names(shp)
+			shp = stars::st_set_dimensions(shp, dimnms[dimid], values = seq_along(geoms))
+			shp = stars::st_set_dimensions(shp, names = dimnms_new)
+		}
+
+		dt = as.data.table(shp)
+
+		if (!is.null(crs) && sf::st_crs(geoms) != crs) {
+			shp = sf::st_transform(shp, crs)
+		} else {
+			shp = geoms
+		}
+		shpTM = shapeTM(shp = shp, tmapID = seq_along(shp), bbox = bbox)
+
+
+		shpclass = "sfc"
+	} else {
+		shp = downsample_stars(shp, max.raster = o$raster.max.cells / (o$fn[1] * o$fn[2]))
+		if (!is.null(crs) && sf::st_crs(shp) != crs) {
+			shp = transwarp(shp, crs, raster.warp = TRUE)
+		}
+
+		dims = stars::st_dimensions(shp)
+		rst = attr(dims, "raster")
+
+		dim_xy = get_xy_dim(shp)
+		dimsxy = dims[names(dim_xy)]
+
+		if (rst$curvilinear) {
+			shp3 = shp
+			attr(shp3, "dimensions")[[rst$dimensions[1]]]$values = 1L:nrow(shp)
+			attr(shp3, "dimensions")[[rst$dimensions[2]]]$values = 1L:ncol(shp)
+			attr(attr(shp3, "dimensions"), "raster")$curvilinear = FALSE
+		} else {
+			shp2 = stars::st_set_dimensions(shp, rst$dimensions[1], values = 1L:nrow(shp))
+			shp3 = stars::st_set_dimensions(shp2, rst$dimensions[2], values = 1L:ncol(shp))
+		}
+
+
+		dt = as.data.table(shp3, center = FALSE)
+
+		if (!all(names(shp3) %in% names(dt))) {
+			subst_names = tail(names(dt), length(shp3))
+			for (i in 1L:length(shp3)) {
+				setnames(dt, subst_names[i], names(shp3)[i])
+			}
+		}
+
+		setnames(dt, names(dim_xy)[1], "X__")
+		setnames(dt, names(dim_xy)[2], "Y__")
+
+		dt[, tmapID__ := as.integer((Y__-1) * nrow(shp) + X__)]
+		dt[, X__:= NULL]
+		dt[, Y__:= NULL]
+
+		data.table::setorder(dt, cols = "tmapID__")
+
+		shp = dimsxy
+		shpclass = "stars"
+
+
+		shpTM = shapeTM(shp = shp, tmapID = 1L:(nrow(shp) * ncol(shp)), bbox = bbox)
+
+	}
+
+	make_by_vars(dt, tmf, smeta)
+
+	#if (is.null(bbox)) bbox = sf::st_bbox(shp)
+
+	dtcols = setdiff(names(dt), "tmapID__")
+
+	filter = if (is.null(filter)) {
+		rep(TRUE, nrow(dt))
+	} else filter[dt$tmapID__]
+	dt[, ':='(sel__ = filter)] # tmapID__ = 1L:nrow(dt),
+
+	for (nm in dtcols) {
+		if (is.factor(dt[[nm]])) {
+			lvls = levels(dt[[nm]])
+			cols = attr(dt[[nm]], "colors")
+			ids = seq_along(lvls)[lvls!="NA"]
+			if (!is.null(cols)) {
+				#if (any(lvls=="NA")) {
+				dt[[nm]] = factor(as.integer(dt[[nm]]), levels = ids, labels = paste(lvls[ids], cols[ids], sep = "=<>="))
+				#}
+				#dt[[nm]] = factor()
+			}
+		}
+	}
+
+	structure(list(shpTM = shpTM, dt = dt, is.main = is.main, dtcols = dtcols, shpclass = shpclass, bbox = bbox, unit = unit, shp_name = shp_name, smeta = smeta), class = "tmapShape")
+}
+
+
+#' @export
+tmapSubsetShp.stars = function(shp, vars) {
+	ids = unique(c(which(names(shp) %in% vars),
+				   which(names(shp) %in% vars)))
+	shp2 = shp[ids]
+	if (!length(vars)) {
+		shp2$dummy__ = TRUE
+	}
+	shp2
+}
+
+
+#' @export
+tmapGetShapeMeta1.stars = function(shp, o) {
+	d = stars::st_dimensions(shp)
+
+	if (!has_raster(shp)) {
+		d_non_xy = local({
+			dimvals = lapply(seq_along(d), function(i) stars::st_get_dimension_values(shp, i))
+			dimsfc = vapply(dimvals, inherits, what = "sfc", FUN.VALUE = logical(1))
+			d[!dimsfc]
+		})
+	} else {
+		d_non_xy = local({
+			dxy = attr(d, "raster")$dimensions
+			d[setdiff(names(d), dxy)]
+		})
+	}
+
+	dims = names(d_non_xy)
+	dims_vals = lapply(dims, function(d) stars::st_get_dimension_values(shp, d))
+	names(dims_vals) = dims
+
+	vars_orig = names(shp)
+	vars = vars_orig #make.names(vars_orig)
+	names(vars) = vars_orig
+	list(vars = vars,
+		 dims = dims,
+		 dims_vals = dims_vals)
+}
+
+#' @export
+tmapGetShapeMeta2.stars = function(shp, smeta, o) {
+	smeta$vars_levs = lapply(seq_len(length(shp)), function(i) {
+		get_fact_levels_na(shp[[i]], o)
+	})
+	names(smeta$vars_levs) = names(shp)
+	smeta
+}
+