From e08e841250edd1c512adbe87baca94a5867a2f06 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Sat, 15 May 2021 08:47:00 -0300 Subject: [PATCH 01/11] add wk implementation --- src/Makevars.in | 1 + src/Makevars.win | 1 + src/wk-impl.c | 2 ++ 3 files changed, 4 insertions(+) create mode 100644 src/wk-impl.c diff --git a/src/Makevars.in b/src/Makevars.in index 11192f02..def4d0c6 100644 --- a/src/Makevars.in +++ b/src/Makevars.in @@ -17,6 +17,7 @@ OBJECTS = cpp-compat.o \ s2-matrix.o \ s2-point.o \ s2-xptr.o \ + wk-impl.o \ s2/base/stringprintf.o \ s2/base/strtoint.o \ s2/encoded_s2cell_id_vector.o \ diff --git a/src/Makevars.win b/src/Makevars.win index c1046ebf..2b532bb8 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -17,6 +17,7 @@ OBJECTS = cpp-compat.o \ s2-matrix.o \ s2-point.o \ s2-xptr.o \ + wk-impl.o \ s2/base/stringprintf.o \ s2/base/strtoint.o \ s2/encoded_s2cell_id_vector.o \ diff --git a/src/wk-impl.c b/src/wk-impl.c new file mode 100644 index 00000000..936224b1 --- /dev/null +++ b/src/wk-impl.c @@ -0,0 +1,2 @@ + +#include "wk-v1-impl.c" From 438439fe185caaa3143df79e0096be9f67fb5c05 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Sun, 16 May 2021 09:42:21 -0300 Subject: [PATCH 02/11] stub tesselating/projecting filter works --- R/wk-utils.R | 20 +++ src/Makevars.in | 2 + src/Makevars.win | 2 + src/RcppExports.cpp | 5 + src/s2-c-api.cpp | 216 ++++++++++++++++++++++ src/wk-c-utils.c | 316 +++++++++++++++++++++++++++++++++ tests/testthat/test-wk-utils.R | 32 ++++ 7 files changed, 593 insertions(+) create mode 100644 R/wk-utils.R create mode 100644 src/s2-c-api.cpp create mode 100644 src/wk-c-utils.c create mode 100644 tests/testthat/test-wk-utils.R diff --git a/R/wk-utils.R b/R/wk-utils.R new file mode 100644 index 00000000..42276bca --- /dev/null +++ b/R/wk-utils.R @@ -0,0 +1,20 @@ + +s2_unprojection_filter <- function(handler, projection = s2_projection_plate_carree(), + tessellate_tol = Inf) { + wk::new_wk_handler( + .Call(c_s2_coord_filter_new, handler, projection, TRUE, tessellate_tol), + subclass = "s2_coord_filter" + ) +} + +s2_projection_filter <- function(handler, projection = s2_projection_plate_carree(), + tessellate_tol = Inf) { + wk::new_wk_handler( + .Call(c_s2_coord_filter_new, handler, projection, FALSE, tessellate_tol), + subclass = "s2_coord_filter" + ) +} + +s2_projection_plate_carree <- function() { + .Call(c_s2_projection_plate_carree) +} diff --git a/src/Makevars.in b/src/Makevars.in index def4d0c6..ff6204df 100644 --- a/src/Makevars.in +++ b/src/Makevars.in @@ -18,6 +18,8 @@ OBJECTS = cpp-compat.o \ s2-point.o \ s2-xptr.o \ wk-impl.o \ + wk-c-utils.o \ + s2-c-api.o \ s2/base/stringprintf.o \ s2/base/strtoint.o \ s2/encoded_s2cell_id_vector.o \ diff --git a/src/Makevars.win b/src/Makevars.win index 2b532bb8..735fb7b3 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -18,6 +18,8 @@ OBJECTS = cpp-compat.o \ s2-point.o \ s2-xptr.o \ wk-impl.o \ + wk-c-utils.o \ + s2-c-api.o \ s2/base/stringprintf.o \ s2/base/strtoint.o \ s2/encoded_s2cell_id_vector.o \ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index da442649..d0e82eac 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1269,6 +1269,9 @@ BEGIN_RCPP END_RCPP } +RcppExport SEXP c_s2_coord_filter_new(SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP c_s2_projection_plate_carree(); + static const R_CallMethodDef CallEntries[] = { {"_s2_cpp_s2_init", (DL_FUNC) &_s2_cpp_s2_init, 0}, {"_s2_cpp_s2_is_collection", (DL_FUNC) &_s2_cpp_s2_is_collection, 1}, @@ -1375,6 +1378,8 @@ static const R_CallMethodDef CallEntries[] = { {"_s2_cpp_s2_buffer_cells", (DL_FUNC) &_s2_cpp_s2_buffer_cells, 4}, {"_s2_s2_xptr_test", (DL_FUNC) &_s2_s2_xptr_test, 1}, {"_s2_s2_xptr_test_op", (DL_FUNC) &_s2_s2_xptr_test_op, 1}, + {"c_s2_coord_filter_new", (DL_FUNC) &c_s2_coord_filter_new, 4}, + {"c_s2_projection_plate_carree", (DL_FUNC) &c_s2_projection_plate_carree, 0}, {NULL, NULL, 0} }; diff --git a/src/s2-c-api.cpp b/src/s2-c-api.cpp new file mode 100644 index 00000000..8ce874cd --- /dev/null +++ b/src/s2-c-api.cpp @@ -0,0 +1,216 @@ + +#include + +#include "s2/s2projections.h" +#include "s2/s2edge_tessellator.h" + +typedef struct { + int (*project)(double* coord_in, double* coord_out, void* data); + int (*unproject)(double* coord_in, double* coord_out, void* data); + double wrap_distance_x; + double wrap_distance_y; + void* data; +} s2_simple_projection_t; + +typedef struct s2_projection_t s2_projection_t; +typedef struct s2_tessellator_t s2_tessellator_t; + +#ifdef __cplusplus +extern "C" { +#endif + +s2_projection_t* s2_projection_create_plate_carree(double scale); +s2_projection_t* s2_projection_create_mercator(double max_x); +s2_projection_t* s2_projection_create(s2_simple_projection_t projection); +void s2_projection_destroy(s2_projection_t* projection); +int s2_projection_project(s2_projection_t* projection, const double* coord_in, double* coord_out); +int s2_projection_unproject(s2_projection_t* projection, const double* coord_in, double* coord_out); + +s2_tessellator_t* s2_tessellator_create(s2_projection_t* projection, double tolerance_radians); +void s2_tessellator_destroy(s2_tessellator_t* tessellator); +int s2_tessellator_reset(s2_tessellator_t* tessellator); +int s2_tessellator_add_r2_point(s2_tessellator_t* tessellator, const double* coord); +int s2_tessellator_add_s2_point(s2_tessellator_t* tessellator, const double* coord); +int s2_tessellator_r2_points_size(s2_tessellator_t* tessellator); +int s2_tessellator_s2_points_size(s2_tessellator_t* tessellator); +int s2_tessellator_r2_point(s2_tessellator_t* tessellator, int i, double* coord); +int s2_tessellator_s2_point(s2_tessellator_t* tessellator, int i, double* coord); + +#ifdef __cplusplus +} +#endif + +// Wrapper class around a projection created using C callables +class SimpleProjection: public S2::Projection { +public: + SimpleProjection(s2_simple_projection_t projection): projection(projection) {} + + R2Point Project(const S2Point& p) const override { + double swap_in[4]; + double swap_out[4]; + swap_in[0] = p.x(); + swap_in[1] = p.y(); + swap_in[2] = p.z(); + projection.project(swap_in, swap_out, projection.data); + return R2Point(swap_out[0], swap_out[1]); + } + + S2Point Unproject(const R2Point& p) const override { + double swap_in[4]; + double swap_out[4]; + swap_in[0] = p.x(); + swap_in[1] = p.y(); + projection.unproject(swap_in, swap_out, projection.data); + return S2Point(swap_out[0], swap_out[1], swap_out[2]); + } + + R2Point FromLatLng(const S2LatLng& ll) const override { + return this->Project(ll.ToPoint()); + } + + S2LatLng ToLatLng(const R2Point& p) const override { + return S2LatLng(this->Unproject(p)); + } + + R2Point wrap_distance() const override { + return R2Point(projection.wrap_distance_x, projection.wrap_distance_y); + } + +private: + s2_simple_projection_t projection; +}; + +s2_projection_t* s2_projection_create(s2_simple_projection_t projection) { + return (s2_projection_t*) new SimpleProjection(projection); +} + +s2_projection_t* s2_projection_create_plate_carree(double scale) { + return (s2_projection_t*) new S2::PlateCarreeProjection(scale); +} + +s2_projection_t* s2_projection_create_mercator(double max_x) { + return (s2_projection_t*) new S2::MercatorProjection(max_x); +} + +void s2_projection_destroy(s2_projection_t* projection) { + if (projection != nullptr) { + delete ((S2::Projection*) projection); + } +} + +int s2_projection_project(s2_projection_t* projection, const double* coord_in, double* coord_out) { + S2Point p(coord_in[0], coord_in[1], coord_in[2]); + R2Point result = ((S2::Projection*) projection)->Project(p); + coord_out[0] = result.x(); + coord_out[1] = result.y(); + return true; +} + +int s2_projection_unproject(s2_projection_t* projection, const double* coord_in, double* coord_out) { + R2Point p(coord_in[0], coord_in[1]); + S2Point result = ((S2::Projection*) projection)->Unproject(p); + coord_out[0] = result.x(); + coord_out[1] = result.y(); + coord_out[2] = result.z(); + return true; +} + +// Wrapper class around a tessellator that also keeps coordinate buffers +// and methods to iterate through them. +class TessellatorWrapper { +public: + TessellatorWrapper(s2_projection_t* projection, double tolerance_radians): + tessellator((S2::Projection*) projection, S1Angle::Radians(tolerance_radians)) {} + + void reset() { + s2points.clear(); + r2points.clear(); + } + + void add_r2_point(const double* coord) { + R2Point pt(coord[0], coord[1]); + if (r2points.size() > 0) { + this->tessellator.AppendUnprojected(r2last, pt, &(this->s2points)); + } + this->r2last = pt; + } + + void add_s2_point(const double* coord) { + S2Point pt(coord[0], coord[1], coord[2]); + if (r2points.size() > 0) { + this->tessellator.AppendProjected(s2last, pt, &(this->r2points)); + } + this->s2last = pt; + } + + int r2_points_size() { + return this->r2points.size(); + } + + void r2_point(int i, double* coord) { + const R2Point& pt = r2points[i]; + coord[0] = pt.x(); + coord[1] = pt.y(); + } + + int s2_points_size() { + return this->s2points.size(); + } + + void s2_point(int i, double* coord) { + const S2Point& pt = s2points[i]; + coord[0] = pt.x(); + coord[1] = pt.y(); + coord[3] = pt.z(); + } + +private: + S2EdgeTessellator tessellator; + std::vector s2points; + std::vector r2points; + R2Point r2last; + S2Point s2last; +}; + +s2_tessellator_t* s2_tessellator_create(s2_projection_t* projection, double tolerance_radians) { + return (s2_tessellator_t*) new TessellatorWrapper(projection, tolerance_radians); +} + +void s2_tessellator_destroy(s2_tessellator_t* tessellator) { + if (tessellator != nullptr) { + delete ((TessellatorWrapper*) tessellator); + } +} + +int s2_tessellator_reset(s2_tessellator_t* tessellator) { + ((TessellatorWrapper*) tessellator)->reset(); + return true; +} + +int s2_tessellator_add_r2_point(s2_tessellator_t* tessellator, const double* coord) { + ((TessellatorWrapper*) tessellator)->add_r2_point(coord); + return true; +} + +int s2_tessellator_add_s2_point(s2_tessellator_t* tessellator, const double* coord) { + ((TessellatorWrapper*) tessellator)->add_s2_point(coord); + return true; +} + +int s2_tessellator_r2_points_size(s2_tessellator_t* tessellator) { + return ((TessellatorWrapper*) tessellator)->r2_points_size(); +} + +int s2_tessellator_s2_points_size(s2_tessellator_t* tessellator) { + return ((TessellatorWrapper*) tessellator)->s2_points_size(); +} + +int s2_tessellator_r2_point(s2_tessellator_t* tessellator, int i, double* coord) { + ((TessellatorWrapper*) tessellator)->r2_point(i, coord); + return true; +} + +int s2_tessellator_s2_point(s2_tessellator_t* tessellator, int i, double* coord) { + ((TessellatorWrapper*) tessellator)->s2_point(i, coord); + return true; +} diff --git a/src/wk-c-utils.c b/src/wk-c-utils.c new file mode 100644 index 00000000..1e46009e --- /dev/null +++ b/src/wk-c-utils.c @@ -0,0 +1,316 @@ + +#define R_NO_REMAP +#include +#include +#include "wk-v1.h" + +typedef struct { + int (*project)(double* coord_in, double* coord_out, void* data); + int (*unproject)(double* coord_in, double* coord_out, void* data); + double wrap_distance_x; + double wrap_distance_y; + void* data; +} s2_simple_projection_t; + +typedef struct s2_projection_t s2_projection_t; +typedef struct s2_tessellator_t s2_tessellator_t; + +#ifdef __cplusplus +extern "C" { +#endif + +s2_projection_t* s2_projection_create_plate_carree(double scale); +s2_projection_t* s2_projection_create_mercator(double max_x); +s2_projection_t* s2_projection_create(s2_simple_projection_t projection); +void s2_projection_destroy(s2_projection_t* projection); +int s2_projection_project(s2_projection_t* projection, const double* coord_in, double* coord_out); +int s2_projection_unproject(s2_projection_t* projection, const double* coord_in, double* coord_out); + +s2_tessellator_t* s2_tessellator_create(s2_projection_t* projection, double tolerance_radians); +void s2_tessellator_destroy(s2_tessellator_t* tessellator); +int s2_tessellator_reset(s2_tessellator_t* tessellator); +int s2_tessellator_add_r2_point(s2_tessellator_t* tessellator, const double* coord); +int s2_tessellator_add_s2_point(s2_tessellator_t* tessellator, const double* coord); +int s2_tessellator_r2_points_size(s2_tessellator_t* tessellator); +int s2_tessellator_s2_points_size(s2_tessellator_t* tessellator); +int s2_tessellator_r2_point(s2_tessellator_t* tessellator, int i, double* coord); +int s2_tessellator_s2_point(s2_tessellator_t* tessellator, int i, double* coord); + +#ifdef __cplusplus +} +#endif + +// Expose projections as external pointers so that they can theoretically be generated +// by other packages. +void s2_projection_xptr_destroy(SEXP projection_xptr) { + s2_projection_destroy((s2_projection_t*) R_ExternalPtrAddr(projection_xptr)); +} + +SEXP c_s2_projection_plate_carree() { + SEXP projection_xptr = PROTECT(R_MakeExternalPtr(s2_projection_create_plate_carree(180), R_NilValue, R_NilValue)); + R_RegisterCFinalizer(projection_xptr, &s2_projection_xptr_destroy); + UNPROTECT(1); + return projection_xptr; +} + +// The s2_coord_filter at the C level is used both for projecting (i.e., +// simple mapping of points from x, y to an x, y, z on the unit sphere) +// and tessellating (i.e., projecting AND adding additional points to +// ensure that "straight" lines in the current projection refer to the same +// line on the surface of the earth when unprojected on to the unit +// sphere. Going the other direction, points are added to ensure that +// great circle edges (S2's assumption) refer to the same line on the +// surface of the earth when exported to the specified projection. +// This is similar to D3's "adaptive resampling" which also handles +// segmentizing and projecting as a single operation. +typedef struct { + s2_projection_t* projection; + s2_tessellator_t* tessellator; + wk_handler_t* next; + wk_meta_t meta_copy; + wk_vector_meta_t vector_meta_copy; + int unproject; + uint32_t coord_id; +} coord_filter_t; + +static inline void modify_meta_for_filter(coord_filter_t* coord_filter, const wk_meta_t* meta) { + memcpy(&(coord_filter->meta_copy), meta, sizeof(wk_meta_t)); + + // if we're projecting we're going to end up with an XY point + // if we're unprojecting, we're going to get an XYZ point + if (coord_filter->unproject) { + coord_filter->meta_copy.flags |= WK_FLAG_HAS_Z; + } else { + coord_filter->meta_copy.flags &= ~WK_FLAG_HAS_Z; + } + + // if we're tesselating a polyline, the 'size' should be set to unknown + if (meta->geometry_type == WK_LINESTRING) { + coord_filter->meta_copy.size = WK_SIZE_UNKNOWN; + } + + // any srid that initially existed is definitely no longer valid + coord_filter->meta_copy.srid = WK_SRID_NONE; +} + +void s2_coord_filter_initialize(int* dirty, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + *dirty = 1; + coord_filter->next->initialize(&coord_filter->next->dirty, coord_filter->next->handler_data); +} + +int s2_coord_filter_vector_start(const wk_vector_meta_t* meta, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + + // if we're projecting we're going to end up with an XY point + // if we're unprojecting, we're going to get an XYZ point + memcpy(&coord_filter->vector_meta_copy, meta, sizeof(wk_vector_meta_t)); + if (coord_filter->unproject) { + coord_filter->vector_meta_copy.flags |= WK_FLAG_HAS_Z; + } else { + coord_filter->vector_meta_copy.flags &= ~WK_FLAG_HAS_Z; + } + + return coord_filter->next->vector_start(&(coord_filter->vector_meta_copy), coord_filter->next->handler_data); +} + +SEXP s2_coord_filter_vector_end(const wk_vector_meta_t* meta, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + // use the modified vector meta from vector_start + return coord_filter->next->vector_end(&(coord_filter->vector_meta_copy), coord_filter->next->handler_data); +} + +int s2_coord_filter_feature_start(const wk_vector_meta_t* meta, R_xlen_t feat_id, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + return coord_filter->next->feature_start(&(coord_filter->vector_meta_copy), feat_id, coord_filter->next->handler_data); +} + +int s2_coord_filter_feature_null(void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + return coord_filter->next->null_feature(coord_filter->next->handler_data); +} + +int s2_coord_filter_feature_end(const wk_vector_meta_t* meta, R_xlen_t feat_id, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + return coord_filter->next->feature_end(&(coord_filter->vector_meta_copy), feat_id, coord_filter->next->handler_data); +} + +int s2_coord_filter_geometry_start(const wk_meta_t* meta, uint32_t part_id, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + modify_meta_for_filter(coord_filter, meta); + if (coord_filter->tessellator) { + s2_tessellator_reset(coord_filter->tessellator); + coord_filter->coord_id = 0; + } + return coord_filter->next->geometry_start(&(coord_filter->meta_copy), part_id, coord_filter->next->handler_data); +} + +int s2_coord_filter_geometry_end(const wk_meta_t* meta, uint32_t part_id, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + modify_meta_for_filter(coord_filter, meta); + return coord_filter->next->geometry_end(&(coord_filter->meta_copy), part_id, coord_filter->next->handler_data); +} + +int s2_coord_filter_ring_start(const wk_meta_t* meta, uint32_t size, uint32_t ring_id, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + modify_meta_for_filter(coord_filter, meta); + if (coord_filter->tessellator) { + s2_tessellator_reset(coord_filter->tessellator); + coord_filter->coord_id = 0; + } + return coord_filter->next->ring_start(&(coord_filter->meta_copy), WK_SIZE_UNKNOWN, ring_id, coord_filter->next->handler_data); +} + +int s2_coord_filter_ring_end(const wk_meta_t* meta, uint32_t size, uint32_t ring_id, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + modify_meta_for_filter(coord_filter, meta); + return coord_filter->next->ring_end(&(coord_filter->meta_copy), WK_SIZE_UNKNOWN, ring_id, coord_filter->next->handler_data); +} + +int s2_coord_filter_error(const char* message, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + return coord_filter->next->error(message, coord_filter->next->handler_data); +} + +void s2_coord_filter_deinitialize(void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + coord_filter->next->deinitialize(coord_filter->next->handler_data); +} + +void s2_coord_filter_finalize(void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + if (coord_filter != NULL) { + s2_tessellator_destroy(coord_filter->tessellator); + // coord_filter->projection is managed using an external pointer + // coord_filter->next is managed using an external pointer + free(coord_filter); + } +} + +// The main show for the coord filter! We define separate functions for the +// projecting and unprojecting case as it's difficult to keep the code +// clean otherwise. +int s2_coord_filter_coord_unproject(const wk_meta_t* meta, const double* coord, uint32_t coord_id, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + modify_meta_for_filter(coord_filter, meta); + double coord_out[4]; + + // if we're not tessellating, skip adding the edge to the tessellator + if ((coord_filter->tessellator == NULL) || (meta->geometry_type == WK_POINT)) { + s2_projection_unproject(coord_filter->projection, coord, coord_out); + return coord_filter->next->coord( + &(coord_filter->meta_copy), + coord_out, coord_id, coord_filter->next->handler_data + ); + } + + // Add this point to the tessellator + s2_tessellator_add_r2_point(coord_filter->tessellator, coord); + + // If we have an edge in the tessellator, reguritate all the points + // except the first one to the next handler. If we have just one point, + // we send it to the handler right away since we'll skip it when + // the next point arrives. + int size = s2_tessellator_s2_points_size(coord_filter->tessellator); + int result; + + if (size == 1) { + s2_projection_unproject(coord_filter->projection, coord, coord_out); + result = coord_filter->next->coord( + &(coord_filter->meta_copy), + coord_out, coord_id, coord_filter->next->handler_data + ); + + if (result != WK_CONTINUE) { + return result; + } + } else if (size >= 2) { + for (int i = 1; i < size; i++) { + s2_tessellator_s2_point(coord_filter->tessellator, i, coord_out); + result = coord_filter->next->coord( + &(coord_filter->meta_copy), + coord_out, coord_filter->coord_id, coord_filter->next->handler_data + ); + coord_filter->coord_id++; + + if (result != WK_CONTINUE) { + return result; + } + } + + // Clear the tessellator and re-add this point to be ready for the next + // point that forms an edge + s2_tessellator_reset(coord_filter->tessellator); + s2_tessellator_add_r2_point(coord_filter->tessellator, coord); + } + + return WK_CONTINUE; +} + +SEXP c_s2_coord_filter_new(SEXP handler_xptr, SEXP projection_xptr, SEXP unproject, SEXP tessellate_tol) { + if (TYPEOF(handler_xptr) != EXTPTRSXP) { + Rf_error("`handler` must be a wk_handler pointer"); + } + + if (TYPEOF(projection_xptr) != EXTPTRSXP) { + Rf_error("`projection` must be an external pointer"); + } + + if (!IS_SIMPLE_SCALAR(unproject, LGLSXP)) { + Rf_error("`unproject` must be TRUE or FALSE"); + } + + if (!IS_SIMPLE_SCALAR(tessellate_tol, REALSXP) || (REAL(tessellate_tol)[0] < 1e-9)) { + Rf_error("`tessellate` must be a double() greter than 1e-9"); + } + + wk_handler_t* handler = wk_handler_create(); + + handler->initialize = &s2_coord_filter_initialize; + handler->vector_start = &s2_coord_filter_vector_start; + handler->vector_end = &s2_coord_filter_vector_end; + + handler->feature_start = &s2_coord_filter_feature_start; + handler->null_feature = &s2_coord_filter_feature_null; + handler->feature_end = &s2_coord_filter_feature_end; + + handler->geometry_start = &s2_coord_filter_geometry_start; + handler->geometry_end = &s2_coord_filter_geometry_end; + + handler->ring_start = &s2_coord_filter_ring_start; + handler->ring_end = &s2_coord_filter_ring_end; + + handler->coord = &s2_coord_filter_coord_unproject; + + handler->error = &s2_coord_filter_error; + + handler->deinitialize = &s2_coord_filter_deinitialize; + handler->finalizer = &s2_coord_filter_finalize; + + coord_filter_t* coord_filter = (coord_filter_t*) malloc(sizeof(coord_filter_t)); + if (coord_filter == NULL) { + wk_handler_destroy(handler); // # nocov + Rf_error("Failed to alloc handler data"); // # nocov + } + + coord_filter->projection = (s2_projection_t*) R_ExternalPtrAddr(projection_xptr); + if (REAL(tessellate_tol)[0] < R_PosInf) { + coord_filter->tessellator = s2_tessellator_create(coord_filter->projection, REAL(tessellate_tol)[0]); + } else { + coord_filter->tessellator = NULL; + } + + coord_filter->unproject = LOGICAL(unproject)[0]; + coord_filter->next = R_ExternalPtrAddr(handler_xptr); + + if (coord_filter->next->api_version != 1) { + Rf_error("Can't run a wk_handler with api_version '%d'", coord_filter->next->api_version); // # nocov + } + + handler->handler_data = coord_filter; + + // include the external pointers as a tag for this external pointer + // which guarnatees that they will not be garbage collected until + // this object is garbage collected + return wk_handler_create_xptr(handler, handler_xptr, projection_xptr); +} diff --git a/tests/testthat/test-wk-utils.R b/tests/testthat/test-wk-utils.R new file mode 100644 index 00000000..df8f60fd --- /dev/null +++ b/tests/testthat/test-wk-utils.R @@ -0,0 +1,32 @@ + +test_that("unprojection using a wk filter works", { + expect_equal( + wk::wk_handle(wk::xy(0, 0), s2_unprojection_filter(wk::xy_writer())), + wk::xyz(1, 0, 0) + ) + + expect_identical( + wk::wk_handle( + wk::wkt("LINESTRING (0 0, 0 45, -60 45)"), + s2_unprojection_filter(wk::wkt_format_handler(precision = 4)) + ), + # 0.7071 ~ sqrt(2) / 2 + "LINESTRING Z (1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071)" + ) + + expect_identical( + wk::wk_handle( + wk::wkt("POLYGON ((0 0, 0 45, -60 45, 0 0))"), + s2_unprojection_filter(wk::wkt_format_handler(precision = 4)) + ), + # 0.7071 ~ sqrt(2) / 2 + "POLYGON Z ((1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071, 1 0 0))" + ) + + expect_identical( + wk::wk_handle(wk::wkt(NA_character_), s2_unprojection_filter(wk::wkt_writer())), + wk::wkt(NA_character_) + ) +}) + + From 569345883e896fe6423c77708966eff3c4325f80 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Sun, 16 May 2021 10:11:46 -0300 Subject: [PATCH 03/11] unprojecting + tessellating --- src/s2-c-api.cpp | 13 ++++++++++--- src/wk-c-utils.c | 9 +++++---- tests/testthat/test-wk-utils.R | 33 +++++++++++++++++++++++++++++++++ 3 files changed, 48 insertions(+), 7 deletions(-) diff --git a/src/s2-c-api.cpp b/src/s2-c-api.cpp index 8ce874cd..ce41250a 100644 --- a/src/s2-c-api.cpp +++ b/src/s2-c-api.cpp @@ -120,27 +120,32 @@ int s2_projection_unproject(s2_projection_t* projection, const double* coord_in, class TessellatorWrapper { public: TessellatorWrapper(s2_projection_t* projection, double tolerance_radians): - tessellator((S2::Projection*) projection, S1Angle::Radians(tolerance_radians)) {} + tessellator((S2::Projection*) projection, S1Angle::Radians(tolerance_radians)), + has_r2_last(false), has_s2_last(false) {} void reset() { s2points.clear(); r2points.clear(); + has_r2_last = false; + has_s2_last = false; } void add_r2_point(const double* coord) { R2Point pt(coord[0], coord[1]); - if (r2points.size() > 0) { + if (has_r2_last) { this->tessellator.AppendUnprojected(r2last, pt, &(this->s2points)); } this->r2last = pt; + this->has_r2_last = true; } void add_s2_point(const double* coord) { S2Point pt(coord[0], coord[1], coord[2]); - if (r2points.size() > 0) { + if (has_s2_last) { this->tessellator.AppendProjected(s2last, pt, &(this->r2points)); } this->s2last = pt; + this->has_s2_last = true; } int r2_points_size() { @@ -170,6 +175,8 @@ class TessellatorWrapper { std::vector r2points; R2Point r2last; S2Point s2last; + bool has_r2_last; + bool has_s2_last; }; s2_tessellator_t* s2_tessellator_create(s2_projection_t* projection, double tolerance_radians) { diff --git a/src/wk-c-utils.c b/src/wk-c-utils.c index 1e46009e..8ea41d9b 100644 --- a/src/wk-c-utils.c +++ b/src/wk-c-utils.c @@ -208,23 +208,24 @@ int s2_coord_filter_coord_unproject(const wk_meta_t* meta, const double* coord, s2_tessellator_add_r2_point(coord_filter->tessellator, coord); // If we have an edge in the tessellator, reguritate all the points - // except the first one to the next handler. If we have just one point, - // we send it to the handler right away since we'll skip it when + // except the first one to the next handler. If we have zero points, + // we send the it to the handler right away since we'll skip it when // the next point arrives. int size = s2_tessellator_s2_points_size(coord_filter->tessellator); int result; - if (size == 1) { + if (size < 2) { s2_projection_unproject(coord_filter->projection, coord, coord_out); result = coord_filter->next->coord( &(coord_filter->meta_copy), coord_out, coord_id, coord_filter->next->handler_data ); + coord_filter->coord_id++; if (result != WK_CONTINUE) { return result; } - } else if (size >= 2) { + } else { for (int i = 1; i < size; i++) { s2_tessellator_s2_point(coord_filter->tessellator, i, coord_out); result = coord_filter->next->coord( diff --git a/tests/testthat/test-wk-utils.R b/tests/testthat/test-wk-utils.R index df8f60fd..bde8e878 100644 --- a/tests/testthat/test-wk-utils.R +++ b/tests/testthat/test-wk-utils.R @@ -29,4 +29,37 @@ test_that("unprojection using a wk filter works", { ) }) +test_that("tessellating + unprojection using a wk filter works", { + # using big examples here, so use a tolerance of 100 km (forces + # adding at least one point) + tol <- 100000 / s2_earth_radius_meters() + expect_equal( + wk::wk_handle( + wk::xy(0, 0), + s2_unprojection_filter(wk::xy_writer(), tessellate_tol = tol) + ), + wk::xyz(1, 0, 0) + ) + + expect_identical( + wk::wk_handle( + wk::wkt("LINESTRING (0 0, 0 45, -60 45)"), + s2_unprojection_filter(wk::wkb_writer(), tessellate_tol = tol) + ) %>% s2_num_points(), + 6L + ) + + expect_identical( + wk::wk_handle( + wk::wkt("POLYGON ((0 0, 0 45, -60 45, 0 0))"), + s2_unprojection_filter(wk::wkb_writer(), tessellate_tol = tol) + ) %>% s2_num_points(), + 8L + ) + + expect_identical( + wk::wk_handle(wk::wkt(NA_character_), s2_unprojection_filter(wk::wkt_writer(), tessellate_tol = tol)), + wk::wkt(NA_character_) + ) +}) From 23b1a64cdecf27987352b64c826bcac711b35f31 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Sun, 16 May 2021 14:02:52 -0300 Subject: [PATCH 04/11] test construction errors --- src/wk-c-utils.c | 2 +- tests/testthat/test-wk-utils.R | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/wk-c-utils.c b/src/wk-c-utils.c index 8ea41d9b..ffedd834 100644 --- a/src/wk-c-utils.c +++ b/src/wk-c-utils.c @@ -258,7 +258,7 @@ SEXP c_s2_coord_filter_new(SEXP handler_xptr, SEXP projection_xptr, SEXP unproje } if (!IS_SIMPLE_SCALAR(unproject, LGLSXP)) { - Rf_error("`unproject` must be TRUE or FALSE"); + Rf_error("`unproject` must be TRUE or FALSE"); // # nocov } if (!IS_SIMPLE_SCALAR(tessellate_tol, REALSXP) || (REAL(tessellate_tol)[0] < 1e-9)) { diff --git a/tests/testthat/test-wk-utils.R b/tests/testthat/test-wk-utils.R index bde8e878..360dc9be 100644 --- a/tests/testthat/test-wk-utils.R +++ b/tests/testthat/test-wk-utils.R @@ -1,4 +1,11 @@ +test_that("the projection/unprojection filter errors for invalid input", { + expect_error(s2_unprojection_filter(NULL), "must be a ") + expect_error(s2_unprojection_filter(wk::wk_void_handler(), projection = NULL), "must be an") + expect_error(s2_unprojection_filter(wk::wk_void_handler(), tessellate_tol = -1), "must be") + expect_error(s2_unprojection_filter(wk::wk_void_handler(), tessellate_tol = -1), "must be") +}) + test_that("unprojection using a wk filter works", { expect_equal( wk::wk_handle(wk::xy(0, 0), s2_unprojection_filter(wk::xy_writer())), From 0c5ecab2605fb741e9df4b3606ae69b941ef59a1 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Sun, 16 May 2021 15:05:41 -0300 Subject: [PATCH 05/11] make sure projection, projection + tessellating works --- src/wk-c-utils.c | 60 ++++++++++++++++++++++++++++++++-- tests/testthat/test-wk-utils.R | 58 ++++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+), 2 deletions(-) diff --git a/src/wk-c-utils.c b/src/wk-c-utils.c index ffedd834..3e10397d 100644 --- a/src/wk-c-utils.c +++ b/src/wk-c-utils.c @@ -248,6 +248,58 @@ int s2_coord_filter_coord_unproject(const wk_meta_t* meta, const double* coord, return WK_CONTINUE; } +int s2_coord_filter_coord_project(const wk_meta_t* meta, const double* coord, uint32_t coord_id, void* handler_data) { + coord_filter_t* coord_filter = (coord_filter_t*) handler_data; + modify_meta_for_filter(coord_filter, meta); + double coord_out[4]; + + if ((coord_filter->tessellator == NULL) || (meta->geometry_type == WK_POINT)) { + s2_projection_project(coord_filter->projection, coord, coord_out); + return coord_filter->next->coord( + &(coord_filter->meta_copy), + coord_out, coord_id, coord_filter->next->handler_data + ); + } + + s2_tessellator_add_s2_point(coord_filter->tessellator, coord); + + int size = s2_tessellator_r2_points_size(coord_filter->tessellator); + int result; + + if (size < 2) { + s2_projection_project(coord_filter->projection, coord, coord_out); + result = coord_filter->next->coord( + &(coord_filter->meta_copy), + coord_out, coord_id, coord_filter->next->handler_data + ); + coord_filter->coord_id++; + + if (result != WK_CONTINUE) { + return result; + } + } else { + for (int i = 1; i < size; i++) { + s2_tessellator_r2_point(coord_filter->tessellator, i, coord_out); + result = coord_filter->next->coord( + &(coord_filter->meta_copy), + coord_out, coord_filter->coord_id, coord_filter->next->handler_data + ); + coord_filter->coord_id++; + + if (result != WK_CONTINUE) { + return result; + } + } + + // Clear the tessellator and re-add this point to be ready for the next + // point that forms an edge + s2_tessellator_reset(coord_filter->tessellator); + s2_tessellator_add_s2_point(coord_filter->tessellator, coord); + } + + return WK_CONTINUE; +} + SEXP c_s2_coord_filter_new(SEXP handler_xptr, SEXP projection_xptr, SEXP unproject, SEXP tessellate_tol) { if (TYPEOF(handler_xptr) != EXTPTRSXP) { Rf_error("`handler` must be a wk_handler pointer"); @@ -281,8 +333,6 @@ SEXP c_s2_coord_filter_new(SEXP handler_xptr, SEXP projection_xptr, SEXP unproje handler->ring_start = &s2_coord_filter_ring_start; handler->ring_end = &s2_coord_filter_ring_end; - handler->coord = &s2_coord_filter_coord_unproject; - handler->error = &s2_coord_filter_error; handler->deinitialize = &s2_coord_filter_deinitialize; @@ -302,6 +352,12 @@ SEXP c_s2_coord_filter_new(SEXP handler_xptr, SEXP projection_xptr, SEXP unproje } coord_filter->unproject = LOGICAL(unproject)[0]; + if (coord_filter->unproject) { + handler->coord = &s2_coord_filter_coord_unproject; + } else { + handler->coord = &s2_coord_filter_coord_project; + } + coord_filter->next = R_ExternalPtrAddr(handler_xptr); if (coord_filter->next->api_version != 1) { diff --git a/tests/testthat/test-wk-utils.R b/tests/testthat/test-wk-utils.R index 360dc9be..62a9267e 100644 --- a/tests/testthat/test-wk-utils.R +++ b/tests/testthat/test-wk-utils.R @@ -70,3 +70,61 @@ test_that("tessellating + unprojection using a wk filter works", { wk::wkt(NA_character_) ) }) + +test_that("projection using a wk filter works", { + expect_equal( + wk::wk_handle(wk::xyz(1, 0, 0), s2_projection_filter(wk::xy_writer())), + wk::xy(0, 0) + ) + + expect_identical( + wk::wk_handle( + wk::wkt("LINESTRING Z (1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071)"), + s2_projection_filter(wk::wkt_format_handler(precision = 4)) + ), + "LINESTRING (0 0, 0 45, -60 45)" + ) + + expect_identical( + wk::wk_handle( + wk::wkt("POLYGON Z ((1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071, 1 0 0))"), + s2_projection_filter(wk::wkt_format_handler(precision = 4)) + ), + "POLYGON ((0 0, 0 45, -60 45, 0 0))" + ) + + expect_identical( + wk::wk_handle(wk::wkt(NA_character_), s2_projection_filter(wk::wkt_writer())), + wk::wkt(NA_character_) + ) +}) + +test_that("projection + tessellating using a wk filter works", { + tol <- 100000 / s2_earth_radius_meters() + + expect_equal( + wk::wk_handle(wk::xyz(1, 0, 0), s2_projection_filter(wk::xy_writer(), tessellate_tol = tol)), + wk::xy(0, 0) + ) + + expect_identical( + wk::wk_handle( + wk::wkt("LINESTRING Z (1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071)"), + s2_projection_filter(wk::wkb_writer(), tessellate_tol = tol) + ) %>% s2_num_points(), + 6L + ) + + expect_identical( + wk::wk_handle( + wk::wkt("POLYGON Z ((1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071, 1 0 0))"), + s2_projection_filter(wk::wkb_writer(), tessellate_tol = tol) + ) %>% s2_num_points(), + 8L + ) + + expect_identical( + wk::wk_handle(wk::wkt(NA_character_), s2_projection_filter(wk::wkt_writer(), tessellate_tol = tol)), + wk::wkt(NA_character_) + ) +}) From bafee855b46a3ba8d8074beeb81147b4973f85c5 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Sun, 16 May 2021 15:16:53 -0300 Subject: [PATCH 06/11] expose mercator projection, drop arbitrary projections for now --- R/wk-utils.R | 4 +++ src/RcppExports.cpp | 2 ++ src/s2-c-api.cpp | 53 ---------------------------------- src/wk-c-utils.c | 16 +++++----- tests/testthat/test-wk-utils.R | 10 +++++++ 5 files changed, 23 insertions(+), 62 deletions(-) diff --git a/R/wk-utils.R b/R/wk-utils.R index 42276bca..f7eaa885 100644 --- a/R/wk-utils.R +++ b/R/wk-utils.R @@ -18,3 +18,7 @@ s2_projection_filter <- function(handler, projection = s2_projection_plate_carre s2_projection_plate_carree <- function() { .Call(c_s2_projection_plate_carree) } + +s2_projection_mercator <- function() { + .Call(c_s2_projection_mercator) +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index d0e82eac..3b1c62e3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1270,6 +1270,7 @@ END_RCPP } RcppExport SEXP c_s2_coord_filter_new(SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP c_s2_projection_mercator(); RcppExport SEXP c_s2_projection_plate_carree(); static const R_CallMethodDef CallEntries[] = { @@ -1379,6 +1380,7 @@ static const R_CallMethodDef CallEntries[] = { {"_s2_s2_xptr_test", (DL_FUNC) &_s2_s2_xptr_test, 1}, {"_s2_s2_xptr_test_op", (DL_FUNC) &_s2_s2_xptr_test_op, 1}, {"c_s2_coord_filter_new", (DL_FUNC) &c_s2_coord_filter_new, 4}, + {"c_s2_projection_mercator", (DL_FUNC) &c_s2_projection_mercator, 0}, {"c_s2_projection_plate_carree", (DL_FUNC) &c_s2_projection_plate_carree, 0}, {NULL, NULL, 0} }; diff --git a/src/s2-c-api.cpp b/src/s2-c-api.cpp index ce41250a..7443e01f 100644 --- a/src/s2-c-api.cpp +++ b/src/s2-c-api.cpp @@ -4,14 +4,6 @@ #include "s2/s2projections.h" #include "s2/s2edge_tessellator.h" -typedef struct { - int (*project)(double* coord_in, double* coord_out, void* data); - int (*unproject)(double* coord_in, double* coord_out, void* data); - double wrap_distance_x; - double wrap_distance_y; - void* data; -} s2_simple_projection_t; - typedef struct s2_projection_t s2_projection_t; typedef struct s2_tessellator_t s2_tessellator_t; @@ -21,7 +13,6 @@ extern "C" { s2_projection_t* s2_projection_create_plate_carree(double scale); s2_projection_t* s2_projection_create_mercator(double max_x); -s2_projection_t* s2_projection_create(s2_simple_projection_t projection); void s2_projection_destroy(s2_projection_t* projection); int s2_projection_project(s2_projection_t* projection, const double* coord_in, double* coord_out); int s2_projection_unproject(s2_projection_t* projection, const double* coord_in, double* coord_out); @@ -40,50 +31,6 @@ int s2_tessellator_s2_point(s2_tessellator_t* tessellator, int i, double* coord) } #endif -// Wrapper class around a projection created using C callables -class SimpleProjection: public S2::Projection { -public: - SimpleProjection(s2_simple_projection_t projection): projection(projection) {} - - R2Point Project(const S2Point& p) const override { - double swap_in[4]; - double swap_out[4]; - swap_in[0] = p.x(); - swap_in[1] = p.y(); - swap_in[2] = p.z(); - projection.project(swap_in, swap_out, projection.data); - return R2Point(swap_out[0], swap_out[1]); - } - - S2Point Unproject(const R2Point& p) const override { - double swap_in[4]; - double swap_out[4]; - swap_in[0] = p.x(); - swap_in[1] = p.y(); - projection.unproject(swap_in, swap_out, projection.data); - return S2Point(swap_out[0], swap_out[1], swap_out[2]); - } - - R2Point FromLatLng(const S2LatLng& ll) const override { - return this->Project(ll.ToPoint()); - } - - S2LatLng ToLatLng(const R2Point& p) const override { - return S2LatLng(this->Unproject(p)); - } - - R2Point wrap_distance() const override { - return R2Point(projection.wrap_distance_x, projection.wrap_distance_y); - } - -private: - s2_simple_projection_t projection; -}; - -s2_projection_t* s2_projection_create(s2_simple_projection_t projection) { - return (s2_projection_t*) new SimpleProjection(projection); -} - s2_projection_t* s2_projection_create_plate_carree(double scale) { return (s2_projection_t*) new S2::PlateCarreeProjection(scale); } diff --git a/src/wk-c-utils.c b/src/wk-c-utils.c index 3e10397d..4c24ca5e 100644 --- a/src/wk-c-utils.c +++ b/src/wk-c-utils.c @@ -4,14 +4,6 @@ #include #include "wk-v1.h" -typedef struct { - int (*project)(double* coord_in, double* coord_out, void* data); - int (*unproject)(double* coord_in, double* coord_out, void* data); - double wrap_distance_x; - double wrap_distance_y; - void* data; -} s2_simple_projection_t; - typedef struct s2_projection_t s2_projection_t; typedef struct s2_tessellator_t s2_tessellator_t; @@ -21,7 +13,6 @@ extern "C" { s2_projection_t* s2_projection_create_plate_carree(double scale); s2_projection_t* s2_projection_create_mercator(double max_x); -s2_projection_t* s2_projection_create(s2_simple_projection_t projection); void s2_projection_destroy(s2_projection_t* projection); int s2_projection_project(s2_projection_t* projection, const double* coord_in, double* coord_out); int s2_projection_unproject(s2_projection_t* projection, const double* coord_in, double* coord_out); @@ -53,6 +44,13 @@ SEXP c_s2_projection_plate_carree() { return projection_xptr; } +SEXP c_s2_projection_mercator() { + SEXP projection_xptr = PROTECT(R_MakeExternalPtr(s2_projection_create_mercator(20037508), R_NilValue, R_NilValue)); + R_RegisterCFinalizer(projection_xptr, &s2_projection_xptr_destroy); + UNPROTECT(1); + return projection_xptr; +} + // The s2_coord_filter at the C level is used both for projecting (i.e., // simple mapping of points from x, y to an x, y, z on the unit sphere) // and tessellating (i.e., projecting AND adding additional points to diff --git a/tests/testthat/test-wk-utils.R b/tests/testthat/test-wk-utils.R index 62a9267e..5724fad2 100644 --- a/tests/testthat/test-wk-utils.R +++ b/tests/testthat/test-wk-utils.R @@ -128,3 +128,13 @@ test_that("projection + tessellating using a wk filter works", { wk::wkt(NA_character_) ) }) + +test_that("mercator projection works", { + expect_equal( + wk::wk_handle( + wk::xy(c(0, 180), 0), + s2_unprojection_filter(s2_projection_filter(wk::xy_writer(), s2_projection_mercator())) + ), + wk::xy(c(0, 20037508), 0) + ) +}) From 021083f802117f35efe648c3bd99001b22baa7b6 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Mon, 17 May 2021 22:38:39 -0300 Subject: [PATCH 07/11] export and document projection/unprojection filters --- NAMESPACE | 4 ++ R/wk-utils.R | 54 ++++++++++++++++++++++++ man/s2_unprojection_filter.Rd | 79 +++++++++++++++++++++++++++++++++++ 3 files changed, 137 insertions(+) create mode 100644 man/s2_unprojection_filter.Rd diff --git a/NAMESPACE b/NAMESPACE index ae8e985c..87ee9ea5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -151,6 +151,9 @@ export(s2_perimeter) export(s2_point) export(s2_project) export(s2_project_normalized) +export(s2_projection_filter) +export(s2_projection_mercator) +export(s2_projection_plate_carree) export(s2_rebuild) export(s2_rebuild_agg) export(s2_simplify) @@ -164,6 +167,7 @@ export(s2_touches) export(s2_touches_matrix) export(s2_union) export(s2_union_agg) +export(s2_unprojection_filter) export(s2_within) export(s2_within_matrix) export(s2_x) diff --git a/R/wk-utils.R b/R/wk-utils.R index f7eaa885..1372a361 100644 --- a/R/wk-utils.R +++ b/R/wk-utils.R @@ -1,4 +1,52 @@ +#' Low-level wk filters and handlers +#' +#' @inheritParams wk::wk_handle +#' @param projection One of [s2_projection_plate_carree()] or +#' [s2_projection_mercator()] +#' @param tessellate_tol An angle in radians. Points will not be added +#' if a line segment is within this distance of a point. +#' +#' @return +#' - `s2_unprojection_filter()`, `s2_projection_filter()`: A `new_wk_handler()` +#' - `s2_projection_plate_carree()`, `s2_projection_mercator()`: An external pointer +#' to an S2 projection. +#' @export +#' +#' @examples +#' library(wk) +#' +#' # simple conversion of individual coordinates *to* unit sphere +#' # space +#' wk_handle( +#' wkt("LINESTRING (0 0, 0 45, -60 45)"), +#' s2_unprojection_filter(wkt_format_handler(5)) +#' ) +#' +#' # simple conversion of individual coordinates *from* unit sphere +#' # space +#' wk_handle( +#' wkt("LINESTRING Z (1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071)"), +#' s2_projection_filter(wkt_format_handler(5)) +#' ) +#' +#' # use tessellate_tol to force points to be added to an edge +#' # unprojection will ensure an edge maintains its cartesian +#' # assumption when transformed to the unit sphere +#' # (i.e., what you probably want when importing a geography) +#' wk_handle( +#' wkt("LINESTRING (0 0, 0 45, -60 45)"), +#' s2_unprojection_filter(wkt_format_handler(5), tessellate_tol = 0.001) +#' ) +#' +#' # projection will ensure an edge maintains its geodesic +#' # assumption when transformed to projected space +#' # (i.e., what you probably want when exporting a geography) +#' wk_handle( +#' wkt("LINESTRING Z (1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071)"), +#' s2_projection_filter(wkt_format_handler(5), tessellate_tol = 0.001) +#' ) +#' s2_unprojection_filter <- function(handler, projection = s2_projection_plate_carree(), tessellate_tol = Inf) { wk::new_wk_handler( @@ -7,6 +55,8 @@ s2_unprojection_filter <- function(handler, projection = s2_projection_plate_car ) } +#' @rdname s2_unprojection_filter +#' @export s2_projection_filter <- function(handler, projection = s2_projection_plate_carree(), tessellate_tol = Inf) { wk::new_wk_handler( @@ -15,10 +65,14 @@ s2_projection_filter <- function(handler, projection = s2_projection_plate_carre ) } +#' @rdname s2_unprojection_filter +#' @export s2_projection_plate_carree <- function() { .Call(c_s2_projection_plate_carree) } +#' @rdname s2_unprojection_filter +#' @export s2_projection_mercator <- function() { .Call(c_s2_projection_mercator) } diff --git a/man/s2_unprojection_filter.Rd b/man/s2_unprojection_filter.Rd new file mode 100644 index 00000000..0844d407 --- /dev/null +++ b/man/s2_unprojection_filter.Rd @@ -0,0 +1,79 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wk-utils.R +\name{s2_unprojection_filter} +\alias{s2_unprojection_filter} +\alias{s2_projection_filter} +\alias{s2_projection_plate_carree} +\alias{s2_projection_mercator} +\title{Low-level wk filters and handlers} +\usage{ +s2_unprojection_filter( + handler, + projection = s2_projection_plate_carree(), + tessellate_tol = Inf +) + +s2_projection_filter( + handler, + projection = s2_projection_plate_carree(), + tessellate_tol = Inf +) + +s2_projection_plate_carree() + +s2_projection_mercator() +} +\arguments{ +\item{handler}{A \link[wk:wk_handle]{wk_handler} object.} + +\item{projection}{One of \code{\link[=s2_projection_plate_carree]{s2_projection_plate_carree()}} or +\code{\link[=s2_projection_mercator]{s2_projection_mercator()}}} + +\item{tessellate_tol}{An angle in radians. Points will not be added +if a line segment is within this distance of a point.} +} +\value{ +\itemize{ +\item \code{s2_unprojection_filter()}, \code{s2_projection_filter()}: A \code{new_wk_handler()} +\item \code{s2_projection_plate_carree()}, \code{s2_projection_mercator()}: An external pointer +to an S2 projection. +} +} +\description{ +Low-level wk filters and handlers +} +\examples{ +library(wk) + +# simple conversion of individual coordinates *to* unit sphere +# space +wk_handle( + wkt("LINESTRING (0 0, 0 45, -60 45)"), + s2_unprojection_filter(wkt_format_handler(5)) +) + +# simple conversion of individual coordinates *from* unit sphere +# space +wk_handle( + wkt("LINESTRING Z (1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071)"), + s2_projection_filter(wkt_format_handler(5)) +) + +# use tessellate_tol to force points to be added to an edge +# unprojection will ensure an edge maintains its cartesian +# assumption when transformed to the unit sphere +# (i.e., what you probably want when importing a geography) +wk_handle( + wkt("LINESTRING (0 0, 0 45, -60 45)"), + s2_unprojection_filter(wkt_format_handler(5), tessellate_tol = 0.001) +) + +# projection will ensure an edge maintains its geodesic +# assumption when transformed to projected space +# (i.e., what you probably want when exporting a geography) +wk_handle( + wkt("LINESTRING Z (1 0 0, 0.7071 0 0.7071, 0.3536 -0.6124 0.7071)"), + s2_projection_filter(wkt_format_handler(5), tessellate_tol = 0.001) +) + +} From da18bde9abaa3b7d62bce40ef81acd738eb14463 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Mon, 17 May 2021 22:46:14 -0300 Subject: [PATCH 08/11] update pkgdown index --- _pkgdown.yml | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index 1c6de209..fb38ec38 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -66,6 +66,10 @@ reference: contents: - s2_closest_feature +- title: Linear Referencing + contents: + - s2_interpolate + - title: Utility Functions contents: - s2_earth_radius_meters @@ -75,3 +79,10 @@ reference: desc: Useful data for testing and demonstrating s2 functions contents: - starts_with("s2_data") + +- title: Low-level Details + desc: Interact with spherical geometry at a low level + contents: + - s2_cell + - s2_cell_is_valid + - s2_unprojection_filter From 7db576a09ce1999193fb9794453cc8f819f2cee3 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Mon, 17 May 2021 23:10:42 -0300 Subject: [PATCH 09/11] fix typo!! --- src/s2-c-api.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/s2-c-api.cpp b/src/s2-c-api.cpp index 7443e01f..28d55e89 100644 --- a/src/s2-c-api.cpp +++ b/src/s2-c-api.cpp @@ -113,7 +113,7 @@ class TessellatorWrapper { const S2Point& pt = s2points[i]; coord[0] = pt.x(); coord[1] = pt.y(); - coord[3] = pt.z(); + coord[2] = pt.z(); } private: From d293918a2c808ccf95f6622cc7fd57faf59cb0f1 Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Thu, 20 May 2021 23:17:55 -0300 Subject: [PATCH 10/11] test early returns for wk filter --- tests/testthat/test-wk-utils.R | 38 ++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tests/testthat/test-wk-utils.R b/tests/testthat/test-wk-utils.R index 5724fad2..98170d6d 100644 --- a/tests/testthat/test-wk-utils.R +++ b/tests/testthat/test-wk-utils.R @@ -129,6 +129,44 @@ test_that("projection + tessellating using a wk filter works", { ) }) +test_that("errors are propagated through the coordinate filter", { + expect_error( + wk::wk_handle(wk::new_wk_wkt("POINT ENTPY"), s2_projection_filter(wk::wk_void_handler())), + "ENTPY" + ) +}) + +test_that("early returns are respected by s2_projection_filter()", { + big_line <- as_wkt(s2::s2_make_line(1:10, 1:10)) + unprojected <- wk::wk_handle(big_line, s2_unprojection_filter(wk::wkt_writer())) + expect_identical( + wk::wk_handle(unprojected, s2_projection_filter(wk::wkt_format_handler(max_coords = 2))), + "LINESTRING (1 1, 2 2...", + ) + expect_identical( + wk::wk_handle( + unprojected, + s2_projection_filter(wk::wkt_format_handler(max_coords = 2), tessellate_tol = 1e-8) + ), + "LINESTRING (1 1, 1.062476 1.062512...", + ) +}) + +test_that("early returns are respected by s2_unprojection_filter()", { + big_line <- as_wkt(s2::s2_make_line(1:10, 1:10)) + expect_identical( + wk::wk_handle(big_line, s2_unprojection_filter(wk::wkt_format_handler(max_coords = 2))), + "LINESTRING Z (0.9996954 0.01744975 0.01745241, 0.998782 0.03487824 0.0348995...", + ) + expect_identical( + wk::wk_handle( + big_line, + s2_unprojection_filter(wk::wkt_format_handler(max_coords = 2), tessellate_tol = 1e-8) + ), + "LINESTRING Z (0.9996954 0.01744975 0.01745241, 0.9996562 0.01853987 0.01854306...", + ) +}) + test_that("mercator projection works", { expect_equal( wk::wk_handle( From c56d3002759caf1a6cd5a62a3117ad3aafe14cce Mon Sep 17 00:00:00 2001 From: Dewey Dunnington Date: Thu, 20 May 2021 23:21:23 -0300 Subject: [PATCH 11/11] NEWS bullet --- NEWS.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/NEWS.md b/NEWS.md index c24cc2ea..76a8b132 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # s2 (development version) +* Added `s2_projection_filter()` and `s2_unprojection_filter()` to + expose the S2 edge tessellator, which can be used to make Cartesian + or great circle assumptions of line segments explicit by adding + points where necessary (#115). * Added an `s2_cell()` vector class to expose a subset of the S2 indexing system to R users (#85, #114). * Added `s2_closest_edges()` to make k-nearest neighbours calculation