diff --git a/R/misc_stars.R b/R/misc_stars.R
index 843d9121..d8fc71af 100644
--- a/R/misc_stars.R
+++ b/R/misc_stars.R
@@ -141,3 +141,13 @@ transwarp = function(x, crs, raster.warp) {
 	if (!is.null(shpcolors)) attr(y[[1]], "colors") = shpcolors
 	y
 }
+
+has_rotate_or_shear = function (x) 
+{
+	dimensions = st_dimensions(x)
+	if (has_raster(x)) {
+		r = attr(dimensions, "raster")
+		!any(is.na(r$affine)) && any(r$affine != 0)
+	}
+	else FALSE
+}
diff --git a/R/step1_helper_facets.R b/R/step1_helper_facets.R
index 8fd5ebc3..bfa8c814 100644
--- a/R/step1_helper_facets.R
+++ b/R/step1_helper_facets.R
@@ -304,28 +304,20 @@ step1_rearrange_facets = function(tmo, o) {
 		
 		smeta = tmapGetShapeMeta2(shp, smeta, c(o, tmg$tmf))
 		if (dev) timing_add(s3 = "get_shape_meta2")
-		
-		
-		
+
 		tmg$tmf = within(tmg$tmf, {
 			
-			
-			
 			gl = list(NULL, NULL, NULL)
 			gn = c(1L, 1L, 1L)
 
-			# assign("gl", gl, envir = .TMAP)
-			# assign("gn", gn, envir = .TMAP)
-			# assign("gisf", is.wrap, envir = .TMAP)
-			
 			for (i in 1L:3L) {
 				byi = get(paste0("by", i))
 				if (!is.null(byi)) {
 					if (byi == "VARS__") {
-						if (!is.null(vl)) gl[[i]] = vl
+						gl[i] = list(vl)
 						gn[i] = vn
 					} else if (byi %in% smeta$vars) {
-						gl[[i]] = smeta$vars_levs[[byi]]
+						gl[i] = list(smeta$vars_levs[[byi]])
 						gn[i] = length(gl[[i]])
 					} else if (byi %in% smeta$dims) {
 						gl[[i]] = smeta$dims_val[[match(byi, smeta$dims)]]	
diff --git a/R/tmapGetShapeMeta.R b/R/tmapGetShapeMeta.R
index b5d29e4e..6c86b757 100644
--- a/R/tmapGetShapeMeta.R
+++ b/R/tmapGetShapeMeta.R
@@ -75,11 +75,15 @@ tmapGetShapeMeta2.stars = function(shp, smeta, o) {
 #' @method tmapGetShapeMeta2 SpatRaster
 #' @export
 tmapGetShapeMeta2.SpatRaster = function(shp, smeta, o) {
-	
-	# slow, needs to be improved with terra functions, e.g. unique and levels
+	if (terra::ncell(shp) > o$raster.max.cells) {
+		# NOTE: this option is not ideal, because categories may be undiscovered
+		# NOTE2: maybe the same should be done with large stars?
+		shp = terra::spatSample(shp, 1e5, method = "regular", as.raster = TRUE)
+	}
 	smeta$vars_levs = lapply(terra::values(shp, dataframe=TRUE), function(dat) {
 		get_fact_levels_na(dat, o)
-	})
+	})		
+
 	names(smeta$vars_levs) = names(shp)
 	smeta
 }
diff --git a/R/tmapShape.R b/R/tmapShape.R
index 9c4acd85..b3fae02b 100644
--- a/R/tmapShape.R
+++ b/R/tmapShape.R
@@ -128,11 +128,14 @@ tmapShape.SpatRaster = function(shp, is.main, crs, bbox, unit, filter, shp_name,
 	
 	#if (is.null(bbox)) bbox = st_bbox(shp)
 	
+	
 	dtcols = setdiff(names(dt), "tmapID__")
 	
 	names(ctabs) = dtcols
 	names(cats) = dtcols
 	
+	make_by_vars(dt, tmf, smeta)
+	
 	if (is.null(filter)) filter = rep(TRUE, nrow(dt))
 	dt[, ':='(sel__ = filter)] # tmapID__ = 1L:nrow(dt), 
 	
diff --git a/sandbox/issues.R b/sandbox/issues.R
index b155eca3..8bf493be 100644
--- a/sandbox/issues.R
+++ b/sandbox/issues.R
@@ -310,6 +310,8 @@ st_crs(s) <- 4326
 st_crs(s0) <- 4326
 st_geotransform(s0) <- c(5, 1.5, 0.2, 0, 0.2, 1.5)
 
+
+
 s0_4326 = st_transform(s0, crs = 4326)
 
 stars:::is_curvilinear(s0)
@@ -449,7 +451,54 @@ tm_shape(nc.32119) + tm_polygons(c("SID74", "SID79"), fill.free = FALSE, fill.le
 
 library(terra)
 library(tmap)
-lc = rast("https://osf.io/download/m92d7")
+lc = rast("sandbox/land_cover.tif") # https://osf.io/download/m92d7")
+
+?terra
+
+levels(lc)
+
+unique(landr)
+
+tm_shape(lc) + tm_raster()
+
+landr = rast(land)
+
+class(landr)
+str(landr)
+terra::lapp(landr, class)
+
+levels(landr)
+
+terra::levels(land)
+
+
+lc2 = terra::spatSample(lc, 1e5, method = "regular", as.raster = T)
+
+tm_shape(lc2) + tm_raster("land_cover") + tm_facets("land_cover")
+tm_shape(lc) + tm_raster("land_cover") + tm_facets("land_cover")
+
+tm_shape(land) + tm_raster("trees") + tm_facets("cover_cls")
+tm_shape(land) + tm_raster("cover_cls") + tm_facets("cover_cls")
+
+
+# to do: 
+
+
+
+
+unique(lc)
+
+x = terra::head(lc, 10)
+
+x = terra::unique(lc)
+
+o = list(drop.NA.facets = TRUE, facet.max = 9)
+
+a = bench::mark(
+	a1 <- get_fact_levels_na(terra::values(lc, dataframe=TRUE)[[1]], o),
+	a2 <- freq(lc),
+	a3 <- unique(lc), check = FALSE
+)
 
 system.time({plot(lc)})