m

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)

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")
 	}
 )

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")
 	}
 )

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
 		}

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}
 }
 

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")

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;
+}
+

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) {