diff --git a/R/extract.R b/R/extract.R
index b50322f4b..3a416806b 100644
--- a/R/extract.R
+++ b/R/extract.R
@@ -47,7 +47,7 @@ getLyrNrs <- function(layer, nms, n) {
 		layer <- match(layer, nms)
 	}
 	if (any(is.na(layer))) {
-		error("extract", "names in argument 'layer' do not match names(x)")
+#		error("extract", "names in argument 'layer' do not match names(x)")
 	}
 	rep_len(layer, n)
 }
@@ -223,6 +223,12 @@ setMethod("extract", signature(x="SpatRaster", y="SpatVector"),
 function(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE, ID=TRUE, weights=FALSE, exact=FALSE, touches=is.lines(y), small=TRUE, layer=NULL, bind=FALSE, raw=FALSE, ...) {
 
 	geo <- geomtype(y)
+	if (!is.null(layer)) {
+		if (length(layer) != nrow(y)) {
+			error("extract", "length(layer) != nrow(y)")
+		}
+	}
+
 	if (geo == "points") {		
 		if (weights || exact) {
 			method <- "bilinear"
@@ -236,12 +242,13 @@ function(x, y, fun=NULL, method="simple", cells=FALSE, xy=FALSE, ID=TRUE, weight
 			if (any(txtfun == "table")) {
 				if (length(fun) > 1) {
 					warn("extract", "'table' cannot be combined with other functions")
-				}			
+				}
 				if (!is.null(layer)) {
 					warn("extract", "argument 'layer' is ignored when 'fun=table'")
 				}
 				e <- extract_table(x, y, ID=ID, weights=weights, exact=exact, touches=touches, small=small, ...)
 			} else {
+				 
 				e <- extract_fun(x, y, txtfun, ID=ID, weights=weights, exact=exact, touches=touches, small=small, bind=bind, layer=layer, ...)
 			}
 			return(e)
diff --git a/R/geom.R b/R/geom.R
index c90beab70..f5f6a53b1 100644
--- a/R/geom.R
+++ b/R/geom.R
@@ -419,8 +419,8 @@ setMethod("spin", signature(x="SpatVector"),
 
 
 setMethod("delaunay", signature(x="SpatVector"),
-	function(x, tolerance=0, as.lines=FALSE) {
-		x@pntr <- x@pntr$delaunay(tolerance, as.lines)
+	function(x, tolerance=0, as.lines=FALSE, constrained=FALSE) {
+		x@pntr <- x@pntr$delaunay(tolerance, as.lines, constrained)
 		messages(x, "delaunay")
 	}
 )
diff --git a/R/merge.R b/R/merge.R
index 7d267a47c..605d32f59 100644
--- a/R/merge.R
+++ b/R/merge.R
@@ -26,20 +26,20 @@ setMethod("merge", signature(x="SpatVector", y="SpatVector"),
 )
 
 setMethod("merge", signature(x="SpatRaster", y="SpatRaster"),
-	function(x, y, ..., first=TRUE, na.rm=TRUE, filename="", overwrite=FALSE, wopt=list()) {
+	function(x, y, ..., first=TRUE, na.rm=TRUE, algo=1, filename="", overwrite=FALSE, wopt=list()) {
 		rc <- sprc(x, y, ...)
 		opt <- spatOptions(filename, overwrite, wopt=wopt)
-		x@pntr <- rc@pntr$merge(first[1], na.rm, opt)
+		x@pntr <- rc@pntr$merge(first[1], na.rm, algo, opt)
 		messages(x, "merge")
 	}
 )
 
 
 setMethod("merge", signature(x="SpatRasterCollection", "missing"),
-	function(x, first=TRUE, na.rm=TRUE, filename="", ...) {
+	function(x, first=TRUE, na.rm=TRUE, algo=1, filename="", ...) {
 		opt <- spatOptions(filename, ...)
 		out <- rast()
-		out@pntr <- x@pntr$merge(first[1], na.rm, opt)
+		out@pntr <- x@pntr$merge(first[1], na.rm, algo, opt)
 		messages(out, "merge")
 	}
 )
diff --git a/R/time.R b/R/time.R
index bd6b08751..e9e0ea76f 100644
--- a/R/time.R
+++ b/R/time.R
@@ -75,6 +75,7 @@ setMethod("time", signature(x="SpatRaster"),
 		d <- x@pntr$time
 		tstep <- x@pntr$timestep
 		if (format != "") {
+			format = match.arg(format, c("seconds", "days", "months", "years", "yearmonths"))
 			if ((format == "months") && (tstep == "years")) {
 				error("time", "cannot extract months from years-time")
 			} else if ((format == "years") && (tstep %in% c("months"))) {
@@ -108,6 +109,9 @@ setMethod("time", signature(x="SpatRaster"),
 		} else if (tstep == "years") {
 			d <- strptime("1970-01-01", "%Y-%m-%d", tz = "UTC") + d
 			as.integer(format(d, "%Y"))
+#		} else if (tstep == "yearweeks") {
+#			d <- strptime("1970-01-01", "%Y-%m-%d", tz = "UTC") + d
+#			yearweek(as.Date(d))
 		} else { # raw
 			d
 		}
diff --git a/man/voronoi.Rd b/man/voronoi.Rd
index 1f4164762..6a82be34f 100644
--- a/man/voronoi.Rd
+++ b/man/voronoi.Rd
@@ -15,7 +15,7 @@ Get a Voronoi diagram or Delaunay triangles for points, or the nodes of lines or
 \usage{
 \S4method{voronoi}{SpatVector}(x, bnd=NULL, tolerance=0, as.lines=FALSE, deldir=FALSE) 
 
-\S4method{delaunay}{SpatVector}(x, tolerance=0, as.lines=FALSE) 
+\S4method{delaunay}{SpatVector}(x, tolerance=0, as.lines=FALSE, constrained=FALSE) 
 }
 
 
@@ -24,6 +24,7 @@ Get a Voronoi diagram or Delaunay triangles for points, or the nodes of lines or
 \item{bnd}{SpatVector to set the outer boundary of the voronoi diagram}
 \item{tolerance}{numeric >= 0, snapping tolerance (0 is no snapping)}
 \item{as.lines}{logical. If \code{TRUE}, lines are returned without the outer boundary}
+\item{constrained}{logical. If \code{TRUE}, a constrained delaunay triangulation is returned}
 \item{deldir}{logical. If \code{TRUE}, the \code{\link[deldir]{deldir}} is used instead of the GEOS C++ library method. It has been reported that \code{deldir} does not choke on very large data sets}
 }
 
diff --git a/src/RcppModule.cpp b/src/RcppModule.cpp
index c94ec61fb..613603c66 100644
--- a/src/RcppModule.cpp
+++ b/src/RcppModule.cpp
@@ -1037,6 +1037,7 @@ RCPP_MODULE(spat){
 		.method("crop", &SpatRasterCollection::crop)
 		.method("addTag", &SpatRasterCollection::addTag)
 		.method("getTags", &SpatRasterCollection::getTags)
+		.method("make_vrt", &SpatRasterCollection::make_vrt)
 	;
 
     class_<SpatRasterStack>("SpatRasterStack")
diff --git a/src/distVector.cpp b/src/distVector.cpp
index 0780624bf..7088f122a 100644
--- a/src/distVector.cpp
+++ b/src/distVector.cpp
@@ -19,6 +19,7 @@
 #include "spatVector.h"
 #include "distance.h"
 #include "geosphere.h"
+#include "geodesic.h"
 #include "vecmath.h"
 
 
@@ -230,7 +231,7 @@ std::vector<double> SpatVector::distance(bool sequential, std::string unit, cons
 		return(d);
 	}
 
-	bool lonlat = is_lonlat(); // m == 0
+	bool lonlat = is_lonlat(); 
 	double m=1;
 	if (!srs.m_dist(m, lonlat, unit)) {
 		setError("invalid unit");
@@ -459,3 +460,151 @@ std::vector<double> SpatVector::pointdistance_seq(const std::vector<double>& px,
 }
 */
 
+
+
+
+void make_dense_lonlat(std::vector<double> &lon, std::vector<double> &lat, const double &interval, const bool &adjust, geod_geodesic &g) {
+	size_t np = lon.size();
+	if (np < 2) {
+		return;
+	}
+	size_t sz = lon.size() * 5;
+	std::vector<double> xout, yout;
+	xout.reserve(sz);
+	yout.reserve(sz);
+	for (size_t i=0; i<(np-1); i++) {
+		if (xout.size() > sz) {
+			sz += (np-i) * 10;
+			xout.reserve(sz);
+			yout.reserve(sz);
+		}
+		double d, azi1, azi2;
+		//double hlat = lat[i] + (lat[i+1] - lat[i])/2;
+		//double hlon = lon[i] + (lon[i+1] - lon[i])/2;
+		//geod_inverse(&g, lat[i], lon[i], hlat, hlon, &d1, &azi1, &azi2);
+		//geod_inverse(&g, hlat, hlon, lat[i+1], lon[i+1], &d2, &azi1, &azi2);
+		//double d = d1 + d2;
+		geod_inverse(&g, lat[i], lon[i], lat[i+1], lon[i+1], &d, &azi1, &azi2);
+		size_t n = floor(d / interval);
+		xout.push_back(lon[i]);
+		yout.push_back(lat[i]);
+		if (n < 2) {
+			continue;
+		}
+		double step = adjust ? d / n : interval;
+		double newlat, newlon;
+		for (size_t j=1; j<n; j++) {
+			geod_direct(&g, lat[i], lon[i], azi1, step*j, &newlat, &newlon, &azi2);
+			// avoid -180 to 180 jumps
+			if ((lon[i] == -180) && (newlon == 180)) newlon = -180;
+			xout.push_back(newlon);
+			yout.push_back(newlat);
+		}
+	}
+	xout.push_back(lon[np-1]);
+	yout.push_back(lat[np-1]);
+	lon = std::move(xout);
+	lat = std::move(yout);
+}
+
+void make_dense_planar(std::vector<double> &x, std::vector<double> &y, double &interval, bool &adjust) {
+	size_t np = x.size();
+	if (np < 2) {
+		return;
+	}
+	size_t sz = x.size() * 5;
+	std::vector<double> xout, yout;
+	xout.reserve(sz);
+	yout.reserve(sz);
+
+	double pi2 = M_PI * 2;
+
+	for (size_t i=0; i<(np-1); i++) {
+		if (xout.size() > sz) {
+			sz += (np-i) * 10;
+			xout.reserve(sz);
+			yout.reserve(sz);
+		}
+		double d = sqrt(pow((x[i+1] - x[i]),2) + pow((y[i+1] - y[i]), 2));
+		size_t n = floor(d / interval);
+		xout.push_back(x[i]);
+		yout.push_back(y[i]);
+		if (n < 2) {
+			continue;
+		}
+
+		double a = fmod(atan2(x[i+1]-x[i], y[i+1]-y[i]), pi2);
+		double step = adjust ? d / n : interval;
+		double distx = step * sin(a);
+		double disty = step * cos(a);
+		for (size_t j=1; j<n; j++) {
+			xout.push_back(x[i] + distx * j);
+			yout.push_back(y[i] + disty * j);
+		}
+	}
+	xout.push_back(x[np-1]);
+	yout.push_back(y[np-1]);
+	x = std::move(xout);
+	y = std::move(yout);
+}
+
+
+
+
+SpatVector SpatVector::densify(double interval, bool adjust, bool ignorelonlat) {
+
+	SpatVector out;
+	if (type() == "points") {
+		out.setError("cannot densify points");
+		return out;
+	}
+	if (interval <= 0) {
+		out.setError("the interval must be > 0");
+		return out;
+	}
+
+	out.srs = srs;
+	if (srs.is_empty()) {
+		out.setError("crs not defined");
+		return(out);
+	}
+	size_t n = size();
+	out.reserve(n);
+	if (is_lonlat() && (!ignorelonlat)) {
+		double a = 6378137.0;
+		double f = 1/298.257223563;
+		struct geod_geodesic geod;
+		geod_init(&geod, a, f);
+
+		for (size_t i=0; i<n; i++) {
+			SpatGeom g = geoms[i];
+			for (size_t j=0; j < geoms[i].size(); j++) {
+				make_dense_lonlat(g.parts[j].x, g.parts[j].y, interval, adjust, geod);
+				if (g.parts[j].hasHoles()) {
+					for (size_t k=0; k < g.parts[j].nHoles(); k++) {
+						make_dense_lonlat(g.parts[j].holes[k].x, g.parts[j].holes[k].y, interval, adjust, geod);
+					}
+				}
+			}
+			g.computeExtent();
+			out.addGeom(g);
+		}
+	} else {
+
+		for (size_t i=0; i<n; i++) {
+			SpatGeom g = geoms[i];
+			for (size_t j=0; j < geoms[i].size(); j++) {
+				make_dense_planar(g.parts[j].x, g.parts[j].y, interval, adjust);
+				if (g.parts[j].hasHoles()) {
+					for (size_t k=0; k < g.parts[j].nHoles(); k++) {
+						make_dense_planar(g.parts[j].holes[k].x, g.parts[j].holes[k].y, interval, adjust);
+					}
+				}
+			}
+			out.addGeom(g);
+		}
+	}
+	out.df = df;
+	return out;
+}
+
diff --git a/src/geos_methods.cpp b/src/geos_methods.cpp
index 0728609a2..c17f51837 100644
--- a/src/geos_methods.cpp
+++ b/src/geos_methods.cpp
@@ -722,6 +722,49 @@ SpatVector SpatVector::shared_paths(SpatVector x, bool index) {
 }
 
 
+bool find_segments(GEOSContextHandle_t hGEOSCtxt, std::vector<double> &x, std::vector<double> &y, std::vector<double> &cx, std::vector<double> &cy, std::vector<size_t> &si, std::vector<double> &sx, std::vector<double> &sy) {
+	size_t n = x.size() - 1;
+	size_t m = cx.size() - 1;
+	double ix, iy;
+	si.resize(0);
+	sx.resize(0);
+	sy.resize(0);
+	for (size_t i=0; i<n; i++) {
+		for (size_t j=0; j<m; j++) {
+			if (GEOSSegmentIntersection_r(hGEOSCtxt, x[i], y[i], x[i+1], cy[j+1], cx[j], cy[j], cx[j+1], cy[j+1], &ix, &iy) == 1) {
+				si.push_back(i);
+				sx.push_back(ix);
+				sy.push_back(iy);
+			}
+		}
+	}
+	return si.size() > 0;
+}
+
+
+SpatVector SpatVector::split_lines(SpatVector v) {
+
+	SpatVector out = *this;
+	std::vector<size_t> si;
+	std::vector<double> sx, sy;
+	GEOSContextHandle_t hGEOSCtxt = geos_init();
+
+	for (size_t i=0; i<v.size(); i++) {
+		SpatVector tmp = v.subset_rows({int(i)});
+		std::vector<int> x = out.relate(tmp, "intersects", true, true);
+		std::vector<std::vector<double>> xy1 = tmp.coordinates();
+		for (size_t i=0; i<x.size(); i++) {
+			if (x[i] == 1) {
+				std::vector<std::vector<double>> xy2 = tmp.coordinates();
+				if (find_segments(hGEOSCtxt, xy1[0], xy1[1], xy2[0], xy2[1], si, sx, sy)) {
+					
+				}
+			}
+		}
+	}
+	return out;
+}	
+
 
 /*
 SpatVector SpatVector::split_polygons(SpatVector lns) {
diff --git a/src/geosphere.cpp b/src/geosphere.cpp
index 671651af8..62ef2e937 100644
--- a/src/geosphere.cpp
+++ b/src/geosphere.cpp
@@ -25,11 +25,11 @@
 #endif
 
 #ifndef M_2PI
-#define M_2PI (3.1415926535897932384626433 * 2.0)
+#define M_2PI (M_PI * 2.0)
 #endif
 
-#ifndef M_PI_2
-#define M_PI_2 (3.1415926535897932384626433 / 2)
+#ifndef M_hPI
+#define M_hPI (M_PI / 2.0)
 #endif
 
 #ifndef WGS84_a
@@ -41,7 +41,7 @@
 #endif
 
 
-#include "Rcpp.h"
+//#include "Rcpp.h"
 
 inline void normLon(double &lon) {
 	lon = fmod(lon + 180, 360.) - 180;
@@ -79,7 +79,7 @@ inline double distance_hav_r(double lon1, double lat1, double lon2, double lat2,
 }
 
 
-
+/*
 double distance_cosdeg(double lon1, double lat1, double lon2, double lat2, double r = 6378137.) {
 	deg2rad(lon1);
 	deg2rad(lon2);
@@ -87,7 +87,7 @@ double distance_cosdeg(double lon1, double lat1, double lon2, double lat2, doubl
 	deg2rad(lat2);
 	return r * acos((sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1-lon2)));
 }
-
+*/
 
 
 void dest_geo(double slon, double slat, double sazi, double dist, double &dlon, double &dlat, double &dazi) {
@@ -287,150 +287,6 @@ std::vector<std::vector<double>> intermediate(double lon1, double lat1, double l
 }
 
 
-void make_dense_lonlat(std::vector<double> &lon, std::vector<double> &lat, const double &interval, const bool &adjust, geod_geodesic &g) {
-	size_t np = lon.size();
-	if (np < 2) {
-		return;
-	}
-	size_t sz = lon.size() * 5;
-	std::vector<double> xout, yout;
-	xout.reserve(sz);
-	yout.reserve(sz);
-	for (size_t i=0; i<(np-1); i++) {
-		if (xout.size() > sz) {
-			sz += (np-i) * 10;
-			xout.reserve(sz);
-			yout.reserve(sz);
-		}
-		double d, azi1, azi2;
-		//double hlat = lat[i] + (lat[i+1] - lat[i])/2;
-		//double hlon = lon[i] + (lon[i+1] - lon[i])/2;
-		//geod_inverse(&g, lat[i], lon[i], hlat, hlon, &d1, &azi1, &azi2);
-		//geod_inverse(&g, hlat, hlon, lat[i+1], lon[i+1], &d2, &azi1, &azi2);
-		//double d = d1 + d2;
-		geod_inverse(&g, lat[i], lon[i], lat[i+1], lon[i+1], &d, &azi1, &azi2);
-		size_t n = floor(d / interval);
-		xout.push_back(lon[i]);
-		yout.push_back(lat[i]);
-		if (n < 2) {
-			continue;
-		}
-		double step = adjust ? d / n : interval;
-		double newlat, newlon;
-		for (size_t j=1; j<n; j++) {
-			geod_direct(&g, lat[i], lon[i], azi1, step*j, &newlat, &newlon, &azi2);
-			// avoid -180 to 180 jumps
-			if ((lon[i] == -180) && (newlon == 180)) newlon = -180;
-			xout.push_back(newlon);
-			yout.push_back(newlat);
-		}
-	}
-	xout.push_back(lon[np-1]);
-	yout.push_back(lat[np-1]);
-	lon = std::move(xout);
-	lat = std::move(yout);
-}
-
-void make_dense_planar(std::vector<double> &x, std::vector<double> &y, double &interval, bool &adjust) {
-	size_t np = x.size();
-	if (np < 2) {
-		return;
-	}
-	size_t sz = x.size() * 5;
-	std::vector<double> xout, yout;
-	xout.reserve(sz);
-	yout.reserve(sz);
-
-	double pi2 = M_PI * 2;
-
-	for (size_t i=0; i<(np-1); i++) {
-		if (xout.size() > sz) {
-			sz += (np-i) * 10;
-			xout.reserve(sz);
-			yout.reserve(sz);
-		}
-		double d = sqrt(pow((x[i+1] - x[i]),2) + pow((y[i+1] - y[i]), 2));
-		size_t n = floor(d / interval);
-		xout.push_back(x[i]);
-		yout.push_back(y[i]);
-		if (n < 2) {
-			continue;
-		}
-
-		double a = fmod(atan2(x[i+1]-x[i], y[i+1]-y[i]), pi2);
-		double step = adjust ? d / n : interval;
-		double distx = step * sin(a);
-		double disty = step * cos(a);
-		for (size_t j=1; j<n; j++) {
-			xout.push_back(x[i] + distx * j);
-			yout.push_back(y[i] + disty * j);
-		}
-	}
-	xout.push_back(x[np-1]);
-	yout.push_back(y[np-1]);
-	x = std::move(xout);
-	y = std::move(yout);
-}
-
-
-SpatVector SpatVector::densify(double interval, bool adjust, bool ignorelonlat) {
-
-	SpatVector out;
-	if (type() == "points") {
-		out.setError("cannot densify points");
-		return out;
-	}
-	if (interval <= 0) {
-		out.setError("the interval must be > 0");
-		return out;
-	}
-
-	out.srs = srs;
-	if (srs.is_empty()) {
-		out.setError("crs not defined");
-		return(out);
-	}
-	size_t n = size();
-	out.reserve(n);
-	if (is_lonlat() && (!ignorelonlat)) {
-		double a = 6378137.0;
-		double f = 1/298.257223563;
-		struct geod_geodesic geod;
-		geod_init(&geod, a, f);
-
-		for (size_t i=0; i<n; i++) {
-			SpatGeom g = geoms[i];
-			for (size_t j=0; j < geoms[i].size(); j++) {
-				make_dense_lonlat(g.parts[j].x, g.parts[j].y, interval, adjust, geod);
-				if (g.parts[j].hasHoles()) {
-					for (size_t k=0; k < g.parts[j].nHoles(); k++) {
-						make_dense_lonlat(g.parts[j].holes[k].x, g.parts[j].holes[k].y, interval, adjust, geod);
-					}
-				}
-			}
-			g.computeExtent();
-			out.addGeom(g);
-		}
-	} else {
-
-		for (size_t i=0; i<n; i++) {
-			SpatGeom g = geoms[i];
-			for (size_t j=0; j < geoms[i].size(); j++) {
-				make_dense_planar(g.parts[j].x, g.parts[j].y, interval, adjust);
-				if (g.parts[j].hasHoles()) {
-					for (size_t k=0; k < g.parts[j].nHoles(); k++) {
-						make_dense_planar(g.parts[j].holes[k].x, g.parts[j].holes[k].y, interval, adjust);
-					}
-				}
-			}
-			out.addGeom(g);
-		}
-	}
-	out.df = df;
-	return out;
-}
-
-
 
 
 std::vector<bool> antipodal(std::vector<double> lon1, std::vector<double> lat1, std::vector<double> lon2, std::vector<double> lat2, double tol=1e-9) {
diff --git a/src/raster_methods.cpp b/src/raster_methods.cpp
index 0decd4f80..163a9153e 100644
--- a/src/raster_methods.cpp
+++ b/src/raster_methods.cpp
@@ -3617,7 +3617,7 @@ bool write_part(SpatRaster& out, SpatRaster& r, const double& hxr, size_t& nl, b
 }
 
 
-SpatRaster SpatRasterCollection::merge(bool first, bool narm, SpatOptions &opt) {
+SpatRaster SpatRasterCollection::merge(bool first, bool narm, int algo, SpatOptions &opt) {
 
 	SpatRaster out;
 	size_t n = size();
@@ -3630,73 +3630,188 @@ SpatRaster SpatRasterCollection::merge(bool first, bool narm, SpatOptions &opt)
 		return(out);
 	}
 
-	std::vector<bool> hvals(n);
-	hvals[0] = ds[0].hasValues();
-	SpatExtent e = ds[0].getExtent();
-	size_t nl = ds[0].nlyr();
-	for (size_t i=1; i<n; i++) {
-		//  lyrs, crs, warncrs, ext, rowcol, res
-		if (!ds[0].compare_geom(ds[i], false, false, opt.get_tolerance(), false, false, false, true)) {
-			out.setError(ds[0].msg.error);
-			return(out);
+	if (algo == 1) {
+		SpatExtent e = ds[0].getExtent();
+		size_t nl = ds[0].nlyr();
+		bool anyvals = false;
+		for (size_t i=1; i<n; i++) {
+			//  lyrs, crs, warncrs, ext, rowcol, res
+			if (!ds[0].compare_geom(ds[i], false, false, opt.get_tolerance(), false, false, false, true)) {
+				out.setError(ds[0].msg.error);
+				return(out);
+			}
+			e.unite(ds[i].getExtent());
+			if (ds[i].hasValues()) anyvals = true;
+			nl = std::max(nl, ds[i].nlyr());
 		}
-		e.unite(ds[i].getExtent());
-		hvals[i] = ds[i].hasValues();
-		nl = std::max(nl, ds[i].nlyr());
-	}
 
-	out = ds[0].geometry(nl, true);
-	//out = ds[0].geometry(1, false);
-	out.setExtent(e, true, true, "");
+		out = ds[0].geometry(nl, true);
+		out.setExtent(e, true, true, "");
+		if (!anyvals) return out;
 
-	bool anyvals = false;
-	for (size_t i=0; i<n; i++) {
-		if (hvals[i]) {
-			anyvals = true;
-			break;
+		double hxr = out.xres()/2;
+
+		std::vector<int> vt = getValueType(true);
+			if (vt.size() == 1) {
+			out.setValueType(vt[0]);
 		}
-	}
-	if (!anyvals) {
-		return out;
-	}
 
-	//out = out.geometry(nl, true);
-	double hxr = out.xres()/2;
+		opt.ncopies = std::max(opt.ncopies, size() + nl);
+		if (!out.writeStart(opt, filenames())) { return out; }
 
-	std::vector<int> vt = getValueType(true);
-		if (vt.size() == 1) {
-		out.setValueType(vt[0]);
-	}
+		std::vector<size_t> seq(n);
+		if (first) {
+			std::iota(seq.rbegin(), seq.rend(), 0);
+		} else {
+			std::iota(seq.begin(), seq.end(), 0);
+		}
 
- 	if (!out.writeStart(opt, filenames())) { return out; }
+		SpatOptions topt(opt);
+		bool warn = false;
+		bool notfirst = false;
+		for (size_t i=0; i<n; i++) {
+			if (!ds[seq[i]].hasValues()) continue;
+			if (narm) {
+				notfirst = i > 0;
+			}
+			if (!write_part(out, ds[seq[i]], hxr, nl, notfirst, warn, topt)) {
+				return out;
+			}
+		}
+		out.writeStop();
+		if (warn) out.addWarning("rasters did not align and were resampled");
 
-	std::vector<size_t> seq(n);
-	if (first) {
-		std::iota(seq.rbegin(), seq.rend(), 0);
-	} else {
-		std::iota(seq.begin(), seq.end(), 0);
-	}
+		return(out);
+		
+	} else if (algo == 2) {
+
+		std::vector<unsigned> use;
+		use.reserve(n);
+		if (ds[0].hasValues()) use.push_back(0);
+		SpatExtent e = ds[0].getExtent();
+		size_t nl = ds[0].nlyr();
+		for (size_t i=1; i<n; i++) {
+			//  lyrs, crs, warncrs, ext, rowcol, res
+			if (!ds[0].compare_geom(ds[i], false, false, opt.get_tolerance(), false, false, false, true)) {
+				out.setError(ds[0].msg.error);
+				return(out);
+			}
+			e.unite(ds[i].getExtent());
+			if (ds[i].hasValues()) use.push_back(i);
+			nl = std::max(nl, ds[i].nlyr());
+		}
 
-	SpatOptions topt(opt);
-	bool warn = false;
-	bool notfirst = false;
-	for (size_t i=0; i<n; i++) {
-		if (!ds[seq[i]].hasValues()) continue;
-		if (narm) {
-			notfirst = i > 0;
+		out = ds[0].geometry(nl, true);
+		out.setExtent(e, true, true, "");
+
+		if (use.empty()) return out;
+
+		std::vector<int> vt = getValueType(true);
+			if (vt.size() == 1) {
+			out.setValueType(vt[0]);
 		}
-		if (!write_part(out, ds[seq[i]], hxr, nl, notfirst, warn, topt)) {
-			return out;
+
+		opt.ncopies = std::max(opt.ncopies, size() + nl);
+		if (!out.writeStart(opt, filenames())) { return out; }
+
+		SpatOptions ropt(opt);
+		
+		for (size_t i=0; i<out.bs.n; i++) {
+			std::vector<std::vector<double>> v;
+			readBlock(out, v, out.bs, i, use, ropt);
+			if (hasError()) {
+				out.writeStop();
+				out.setError(getError());
+				return out;
+			}
+			std::vector<size_t> sizes(v.size());
+
+			size_t n = nl * out.bs.nrows[i] * out.ncol();
+			size_t m = v.size();
+
+			bool msz = false;
+			std::vector<size_t> sz(m);
+			for (size_t j=0; j<m; j++) {
+				sz[j] = v[j].size();
+				if (v[j].size() < n) msz = true;
+				if (v[j].size() > n) {
+				out.setError("something is not right. Exected: " + std::to_string(n) + " got: " + std::to_string(v[j].size()) + " values");
+					return out;
+				}
+			}
+
+			if (m == 1) {
+				if (!out.writeBlock(v[0], i)) return out;				
+			} else if (first) {
+				recycle(v[0], n);
+				if (msz) { // with recycling
+					for (size_t j=0; j<n; j++) {
+						for (size_t k=1; k<m; k++) {
+							if (std::isnan(v[0][j])) {
+								v[0][j] = v[k][j%sz[k]];
+							} else {
+								continue;
+							}
+						}
+					}
+				} else {
+					for (size_t j=0; j<n; j++) {
+						for (size_t k=1; k<m; k++) {
+							if (std::isnan(v[0][j])) {
+								v[0][j] = v[k][j];
+							} else {
+								continue;
+							}
+						}
+					}
+				}
+				if (!out.writeBlock(v[0], i)) return out;	
+			} else { // last
+				m -= 1;
+				recycle(v[m], n);
+				if (msz) {
+					for (size_t j=0; j<n; j++) {
+						for (long k=(m-1); k >= 0; k--) {
+							if (std::isnan(v[m][j])) {
+								v[m][j] = v[k][j%sz[k]];
+							} else {
+								continue;
+							}
+						}
+					}
+				} else {
+					for (size_t j=0; j<n; j++) {
+						for (long k=(m-1); k >= 0; k--) {
+							if (std::isnan(v[m][j])) {
+								v[m][j] = v[k][j];
+							} else {
+								continue;
+							}
+						}
+					}
+				}
+				if (!out.writeBlock(v[m], i)) return out;
+			}
 		}
+		out.writeStop();
+		return(out);
+		
+	} else if (algo==3) {
+		
+		SpatOptions vopt(opt);
+		std::string fname = make_vrt({}, !first, vopt);
+		SpatRaster v(fname, {}, {}, {}, {});
+		return v.writeRaster(opt);
+		
+	} else {
+		out.setError("invalid algo (should be 1, 2, or 3)");
+		return(out);		
 	}
-	out.writeStop();
-	if (warn) out.addWarning("rasters did not align and were resampled");
-
-	return(out);
 }
 
 
 
+
 bool overlaps(const std::vector<unsigned>& r1, const std::vector<unsigned>& r2, 
 			  const std::vector<unsigned>& c1, const std::vector<unsigned>& c2) {
 	size_t n = r1.size();
@@ -3720,12 +3835,10 @@ SpatRaster SpatRasterCollection::mosaic(std::string fun, SpatOptions &opt) {
 		return out;
 	}
 	
-	if (fun == "first") {
-		return merge(true, true, opt);
-	}
-	if (fun == "last") {
-		return merge(false, true, opt);
+	if ((fun == "first") || (fun == "last")) {
+		return merge(fun=="first", true, 1, opt);
 	}
+
 	size_t n = size();
 
 	if (n == 0) {
diff --git a/src/spatRasterMultiple.cpp b/src/spatRasterMultiple.cpp
index 6f80ec955..27a172689 100644
--- a/src/spatRasterMultiple.cpp
+++ b/src/spatRasterMultiple.cpp
@@ -17,6 +17,7 @@
 
 #include "spatRasterMultiple.h"
 #include "string_utils.h"
+#include "file_utils.h"
 
 
 SpatRasterCollection SpatRasterCollection::deepCopy() { return *this; }
@@ -97,6 +98,55 @@ void SpatRasterCollection::erase(size_t i) {
 	}
 }
 
+
+std::string SpatRasterCollection::make_vrt(std::vector<std::string> options, bool reverse, SpatOptions &opt) {
+
+	std::vector<std::string> ff = filenames();
+	SpatOptions xopt(opt);
+	for (size_t i=0; i<ff.size(); i++) {
+		if (ff[i] == "") {
+			ff[i] = tempFile(xopt.get_tempdir(), xopt.tmpfile, "_temp_raster.tif");
+			xopt.set_filenames({ff[i]});
+			SpatRaster out = ds[i].writeRaster(xopt);
+			if (out.hasError()) {
+				setError(out.getError());
+				return "";
+			}
+		}
+	}
+	SpatRaster tmp;
+	if (reverse) std::reverse(ff.begin(), ff.end());
+	return tmp.make_vrt(ff, options, opt);
+}
+
+
+void SpatRasterCollection::readBlock(SpatRaster &r, std::vector<std::vector<double>> &v, BlockSize bs, size_t i, std::vector<unsigned> use, SpatOptions opt){
+
+	if ((bs.row[i] + bs.nrows[i]) > r.nrow()) {
+		setError("invalid rows/columns");
+		return;
+	}
+	if (bs.nrows[i]==0) {
+		return;
+	}
+	SpatExtent re = r.getExtent();
+	double yres = r.yres();
+	double ymx = re.ymax - bs.row[i] * yres;
+	double ymn = re.ymax - (bs.row[i] + bs.nrows[i]) * yres;
+	SpatExtent e = {re.xmin, re.xmax, ymn, ymx};
+	SpatRasterCollection x = crop(e, "near", true, use, opt);
+	if (x.hasError()) {
+		setError(x.getError());
+		return;
+	}
+	v.resize(x.size());
+	for (size_t i=0; i< x.size(); i++) {
+		x.ds[i].readValues(v[i], 0, x.ds[i].nrow(), 0, x.ds[i].ncol());
+	}
+}
+
+
+
 SpatRasterCollection SpatRasterCollection::crop(SpatExtent e, std::string snap, bool expand, std::vector<unsigned> use, SpatOptions &opt) {
 
 	SpatRasterCollection out;
@@ -114,9 +164,11 @@ SpatRasterCollection SpatRasterCollection::crop(SpatExtent e, std::string snap,
 			SpatExtent xe = e.intersect(ds[i].getExtent());
 			if (xe.valid_notempty()) {
 				SpatRaster r = ds[i].crop(e, snap, expand, ops);
-				if (!r.hasError()) {
-					out.push_back(r, "");
-				}
+				if (r.hasError()) {
+					out.setError(r.getError());
+					return out;
+				} 
+				out.push_back(r, "");
 			}
 		}
 	} else {
@@ -124,9 +176,11 @@ SpatRasterCollection SpatRasterCollection::crop(SpatExtent e, std::string snap,
 			SpatExtent xe = e.intersect(ds[use[i]].getExtent());
 			if (xe.valid_notempty()) {
 				SpatRaster r = ds[use[i]].crop(e, snap, expand, ops);
-				if (!r.hasError()) {
-					out.push_back(r, "");
+				if (r.hasError()) {
+					out.setError(r.getError());
+					return out;
 				}
+				out.push_back(r, "");
 			}
 		}
 	}
diff --git a/src/spatRasterMultiple.h b/src/spatRasterMultiple.h
index 825f73d17..f5c36ab06 100644
--- a/src/spatRasterMultiple.h
+++ b/src/spatRasterMultiple.h
@@ -125,11 +125,15 @@ class SpatRasterCollection {
 		void resize(size_t n);
 		void push_back(SpatRaster r, std::string name);
 		void erase(size_t i); 
+		
+		void readBlock(SpatRaster &r, std::vector<std::vector<double>> &v, BlockSize bs, size_t i, std::vector<unsigned> use, SpatOptions opt);
+		std::string make_vrt(std::vector<std::string> options, bool reverse, SpatOptions &opt);
+		
 		SpatRasterCollection crop(SpatExtent e, std::string snap, bool expand, std::vector<unsigned> use, SpatOptions &opt);
 		SpatRasterCollection cropmask(SpatVector v, std::string snap, bool touches, bool expand, std::vector<unsigned> use, SpatOptions &opt);
 		std::vector<int> getValueType(bool unique);
 
-		SpatRaster merge(bool first, bool narm, SpatOptions &opt);
+		SpatRaster merge(bool first, bool narm, int algo, SpatOptions &opt);
 		SpatRaster morph(SpatRaster &x, SpatOptions &opt);
 		SpatRaster mosaic(std::string fun, SpatOptions &opt);
 		SpatRaster summary(std::string fun, SpatOptions &opt);
diff --git a/src/spatVector.h b/src/spatVector.h
index 7afcfa36f..c4605b496 100644
--- a/src/spatVector.h
+++ b/src/spatVector.h
@@ -328,7 +328,7 @@ class SpatVector {
 		SpatVector snap(double tolerance);
 		SpatVector snapto(SpatVector y, double tolerance);
 		SpatVector thin(double threshold);
-
+		SpatVector split_lines(SpatVector v);
 		SpatVector allerretour();
 		SpatVectorCollection bienvenue();
 		SpatVector aggregate(bool dissolve);