x

rspatial · Apr 18, 2021 · 622a6e2 · 622a6e2
1 parent 7d1cf85
commit 622a6e2
Show file tree

Hide file tree

Showing 7 changed files with 141 additions and 138 deletions.
diff --git a/R/extractPoints.R b/R/extractPoints.R
@@ -50,6 +50,7 @@ function(x, y, ..., df=FALSE, sp=FALSE){
 
 .xyValues <- function(object, xy, method='simple', buffer=NULL, small=FALSE, cellnumbers=FALSE, fun=NULL, na.rm=TRUE, layer, nl, df=FALSE, factors=FALSE, sp=FALSE, ...) { 
 
+
 	nlyrs <- nlayers(object)
 	if (nlyrs > 1) {
 		if (missing(layer)) { layer <- 1 } 

diff --git a/R/xyValuesBuffer.R b/R/xyValuesBuffer.R
@@ -200,6 +200,8 @@
 	}
 
 	if (! is.null(fun)) {
+
+		fun <- match.fun(fun)
 		if (na.rm) {
 			fun2 <- function(x){
 						x <- stats::na.omit(x)

diff --git a/src/area.cpp b/src/area.cpp
diff --git a/src/area.h b/src/area.h
diff --git a/src/distance.cpp b/src/distance.cpp
@@ -276,3 +276,132 @@ std::vector<std::vector<double> > destpoint_plane(std::vector<double>  x, std::v
 	return(out);
 }
 
+
+
+
+double area_polygon_lonlat(std::vector<double> lon, std::vector<double> lat, double a, double f) {
+	struct geod_geodesic g;
+	struct geod_polygon p;
+	geod_init(&g, a, f);
+	geod_polygon_init(&p, 0);
+	int n = lat.size();
+	for (int i=0; i < n; i++) {
+		geod_polygon_addpoint(&g, &p, lat[i], lon[i]);
+	}
+	double area, P;
+	geod_polygon_compute(&g, &p, 0, 1, &area, &P);
+	return(area < 0 ? -area : area);
+}
+
+std::vector<double> area_polygon_lonlat(std::vector<double> lon, std::vector<double> lat, std::vector<int> pols, std::vector<int> parts, std::vector<int> holes, double a, double f) {
+
+	std::vector<double> out;
+	struct geod_geodesic g;
+	struct geod_polygon p;
+	geod_init(&g, a, f);
+	geod_polygon_init(&p, 0);
+
+	double area, P, pa, tota;
+	int pol = 1;
+	int part = 1;
+	int n = lon.size();
+	tota = 0;
+	for (int i=0; i < n; i++) {
+		if (parts[i] != part || pols[i] != pol) {
+			geod_polygon_compute(&g, &p, 0, 1, &area, &P);
+			pa = fabs(area);
+			tota += (holes[i-1] > 0 ? -pa : pa); // hole
+			part = parts[i]; 
+			if (pols[i] != pol) {
+				out.push_back(tota);
+				tota = 0;
+				pol = pols[i];
+			}
+			geod_polygon_init(&p, 0);
+		} 
+		geod_polygon_addpoint(&g, &p, lat[i], lon[i]);		
+	}
+	geod_polygon_compute(&g, &p, 0, 1, &area, &P);
+	pa = fabs(area);
+	tota += (holes[n-1] > 0 ? -pa : pa); // hole
+	out.push_back(tota);
+	return(out);
+}
+
+
+
+
+double area_polygon_plane(std::vector<double> x, std::vector<double> y) {
+// based on http://paulbourke.net/geometry/polygonmesh/source1.c
+	int n = x.size();
+	double area = x[n-1] * y[0];
+	area -= y[n-1] * x[0];
+	for (int i=0; i < (n-1); i++) {
+		area += x[i] * y[i+1];
+		area -= x[i+1] * y[i];
+	}
+	area /= 2;
+	return(area < 0 ? -area : area);
+}
+
+/*
+std::vector<double> area_polygon_plane(std::vector<double> x, std::vector<double> y, std::vector<int> pols, std::vector<int> parts, std::vector<int> holes) {
+
+	std::vector<double> out;
+	std::vector<double> px;
+	std::vector<double> py;
+	int pol = 1;
+	int part = 1;
+	int n = x.size();
+	double tota = 0;
+	double pa;
+	for (int i=0; i < n; i++) {
+		if (parts[i] != part || pols[i] != pol) {
+			pa = area_polygon_plane(px, py);
+			tota += (holes[i-1] > 0 ? -pa : pa);
+			part = parts[i];
+			if (pols[i] != pol) {
+				out.push_back(tota);
+				tota = 0;
+				pol = pols[i];
+			}
+			px.resize(0);
+			py.resize(0);
+		} 
+		px.push_back(x[i]);	
+		py.push_back(y[i]);				
+	}
+	pa = area_polygon_plane(px, py);
+	tota += (holes[n-1] > 0 ? -pa : pa);
+	out.push_back(tota);
+	return(out);
+}
+*/
+
+std::vector<double> area_polygon_plane(std::vector<double> x, std::vector<double> y, std::vector<int> pols, std::vector<int> parts, std::vector<int> holes) {
+
+	std::vector<double> out;
+	int pol = 1;
+	int part = 1;
+	int n = x.size();
+	double tota = 0;
+	double pa;
+	int ps = 0;
+	for (int i=0; i < n; i++) {
+		if (parts[i] != part || pols[i] != pol) {
+			pa = area_polygon_plane(std::vector<double> (x.begin() + ps, x.begin() + i - 1), std::vector<double> (y.begin() + ps, y.begin() + i - 1));
+			tota += (holes[i-1] > 0 ? -pa : pa);
+			part = parts[i];
+			ps = i;
+			if (pols[i] != pol) {
+				out.push_back(tota);
+				tota = 0;
+				pol = pols[i];
+			}
+		} 
+	}
+	pa = area_polygon_plane(std::vector<double> (x.begin() + ps, x.end()), std::vector<double> (y.begin() + ps, y.end()));
+	tota += (holes[n-1] > 0 ? -pa : pa);
+	out.push_back(tota);
+	return(out);
+}
diff --git a/src/distance.h b/src/distance.h
@@ -22,3 +22,12 @@ std::vector<std::vector<double> > destpoint_lonlat(std::vector<double> longitude
 
 std::vector<double> destpoint_plane(double x, double y, double bearing, double distance);
 std::vector<std::vector<double> > destpoint_plane(std::vector<double>  x, std::vector<double>  y, std::vector<double>  bearing, std::vector<double>  distance);
+
+
+// area 
+double area_polygon_lonlat(std::vector<double> lon, std::vector<double> lat, double a, double f);
+double area_polygon_plane(std::vector<double> x, std::vector<double> y);
+
+std::vector<double> area_polygon_lonlat(std::vector<double> lon, std::vector<double> lat, std::vector<int> pols, std::vector<int> parts, std::vector<int> holes, double a, double f);
+std::vector<double> area_polygon_plane(std::vector<double> x, std::vector<double> y, std::vector<int> pols, std::vector<int> parts, std::vector<int> holes);
+
diff --git a/src/raster_distance.cpp b/src/raster_distance.cpp
@@ -5,7 +5,6 @@ July 2016
 #include <Rcpp.h>
 
 #include "distance.h"
-#include "area.h"
 
 
 // [[Rcpp::export(name = ".get_area_polygon")]]