diff --git a/NAMESPACE b/NAMESPACE index 0e8d30da..02390de1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -73,6 +73,7 @@ S3method(str,s2_cell_union) S3method(str,s2_xptr) S3method(unique,s2_cell) S3method(unlist,s2_cell_union) +S3method(wk_handle,s2_geography) export(as_s2_cell) export(as_s2_cell_union) export(as_s2_geography) @@ -202,4 +203,5 @@ importFrom(Rcpp,sourceCpp) importFrom(utils,str) importFrom(wk,as_wkb) importFrom(wk,as_wkt) +importFrom(wk,wk_handle) useDynLib(s2, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index 7ed3ac9b..06d3a941 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -253,42 +253,10 @@ cpp_s2_cell_common_ancestor_level_agg <- function(cellId) { .Call(`_s2_cpp_s2_cell_common_ancestor_level_agg`, cellId) } -cpp_s2_geog_point <- function(x, y) { - .Call(`_s2_cpp_s2_geog_point`, x, y) -} - -cpp_s2_make_line <- function(x, y, featureId) { - .Call(`_s2_cpp_s2_make_line`, x, y, featureId) -} - -cpp_s2_make_polygon <- function(x, y, featureId, ringId, oriented, check) { - .Call(`_s2_cpp_s2_make_polygon`, x, y, featureId, ringId, oriented, check) -} - -s2_geography_from_wkb <- function(wkb, oriented, check) { - .Call(`_s2_s2_geography_from_wkb`, wkb, oriented, check) -} - -s2_geography_from_wkt <- function(wkt, oriented, check) { - .Call(`_s2_s2_geography_from_wkt`, wkt, oriented, check) -} - s2_geography_full <- function(x) { .Call(`_s2_s2_geography_full`, x) } -s2_geography_to_wkt <- function(s2_geography, precision, trim) { - .Call(`_s2_s2_geography_to_wkt`, s2_geography, precision, trim) -} - -s2_geography_to_wkb <- function(s2_geography, endian) { - .Call(`_s2_s2_geography_to_wkb`, s2_geography, endian) -} - -s2_geography_format <- function(s2_geography, maxCoords, precision, trim) { - .Call(`_s2_s2_geography_format`, s2_geography, maxCoords, precision, trim) -} - s2_lnglat_from_numeric <- function(lng, lat) { .Call(`_s2_s2_lnglat_from_numeric`, lng, lat) } diff --git a/R/data.R b/R/data.R index 9bd06bb1..23ee0ea1 100644 --- a/R/data.R +++ b/R/data.R @@ -70,3 +70,10 @@ s2_data_cities <- function(name = NULL) { as_s2_geography(wkb) } + +#' Example Geometries +#' +#' These geometries are toy examples useful for testing various coordinate +#' shuffling operations in the s2 package. +#' +"s2_data_example_wkt" diff --git a/R/s2-constructors-formatters.R b/R/s2-constructors-formatters.R index 8926f67a..569f13c5 100644 --- a/R/s2-constructors-formatters.R +++ b/R/s2-constructors-formatters.R @@ -70,31 +70,32 @@ #' s2_as_binary(geog) #' s2_geog_point <- function(longitude, latitude) { - recycled <- recycle_common(longitude, latitude) - new_s2_xptr(cpp_s2_geog_point(recycled[[1]], recycled[[2]]), "s2_geography") + wk::wk_handle(wk::xy(longitude, latitude), s2_geography_writer()) } #' @rdname s2_geog_point #' @export s2_make_line <- function(longitude, latitude, feature_id = 1L) { - recycled <- recycle_common(longitude, latitude, feature_id) - new_s2_xptr(cpp_s2_make_line(recycled[[1]], recycled[[2]], featureId = recycled[[3]]), "s2_geography") + wk::wk_handle( + wk::xy(longitude, latitude), + wk::wk_linestring_filter( + s2_geography_writer(), + feature_id = as.integer(feature_id) + ) + ) } #' @rdname s2_geog_point #' @export s2_make_polygon <- function(longitude, latitude, feature_id = 1L, ring_id = 1L, oriented = FALSE, check = TRUE) { - recycled <- recycle_common(longitude, latitude, feature_id, ring_id) - new_s2_xptr( - cpp_s2_make_polygon( - recycled[[1]], recycled[[2]], - featureId = recycled[[3]], - ringId = recycled[[4]], - oriented = oriented, - check = check - ), - "s2_geography" + wk::wk_handle( + wk::xy(longitude, latitude), + wk::wk_polygon_filter( + s2_geography_writer(oriented = oriented, check = check), + feature_id = as.integer(feature_id), + ring_id = as.integer(ring_id) + ) ) } @@ -102,14 +103,12 @@ s2_make_polygon <- function(longitude, latitude, feature_id = 1L, ring_id = 1L, #' @export s2_geog_from_text <- function(wkt_string, oriented = FALSE, check = TRUE) { attributes(wkt_string) <- NULL - wk::validate_wk_wkt(wk::new_wk_wkt(wkt_string)) - new_s2_xptr( - s2_geography_from_wkt( - wkt_string, - oriented = oriented, - check = check - ), - "s2_geography" + wkt <- wk::new_wk_wkt(wkt_string, geodesic = TRUE) + wk::validate_wk_wkt(wkt) + + wk::wk_handle( + wkt, + s2_geography_writer(oriented = oriented, check = check) ) } @@ -117,25 +116,31 @@ s2_geog_from_text <- function(wkt_string, oriented = FALSE, check = TRUE) { #' @export s2_geog_from_wkb <- function(wkb_bytes, oriented = FALSE, check = TRUE) { attributes(wkb_bytes) <- NULL - wk::validate_wk_wkb(wk::new_wk_wkb(wkb_bytes)) - new_s2_xptr( - s2_geography_from_wkb( - wkb_bytes, - oriented = oriented, - check = check - ), - "s2_geography" + wkb <- wk::new_wk_wkb(wkb_bytes) + wk::validate_wk_wkb(wkb) + wk::wk_handle( + wkb, + s2_geography_writer(oriented = oriented, check = check) ) } #' @rdname s2_geog_point #' @export s2_as_text <- function(x, precision = 16, trim = TRUE) { - s2_geography_to_wkt(as_s2_geography(x), precision = precision, trim = trim) + wkt <- wk::wk_handle( + as_s2_geography(x), + wk::wkt_writer(precision = precision, trim = trim) + ) + + attributes(wkt) <- NULL + wkt } #' @rdname s2_geog_point #' @export s2_as_binary <- function(x, endian = wk::wk_platform_endian()) { - structure(s2_geography_to_wkb(as_s2_geography(x), endian = endian), class = "blob") + structure( + wk::wk_handle(as_s2_geography(x), wk::wkb_writer(endian = endian)), + class = "blob" + ) } diff --git a/R/s2-geography.R b/R/s2-geography.R index c89684c3..5a2280ba 100644 --- a/R/s2-geography.R +++ b/R/s2-geography.R @@ -48,7 +48,7 @@ as_s2_geography.s2_geography <- function(x, ...) { #' @export as_s2_geography.s2_lnglat <- function(x, ...) { df <- data_frame_from_s2_lnglat(x) - new_s2_xptr(cpp_s2_geog_point(df[[1]], df[[2]]), "s2_geography") + s2_geog_point(df[[1]], df[[2]]) } #' @rdname as_s2_geography @@ -76,28 +76,22 @@ as_s2_geography.wk_wkb <- function(x, ..., oriented = FALSE, check = TRUE) { } } - new_s2_xptr( - s2_geography_from_wkb(x, oriented = oriented, check = check), - "s2_geography" + wk::wk_handle( + x, + s2_geography_writer(oriented = oriented, check = check) ) } #' @rdname as_s2_geography #' @export as_s2_geography.WKB <- function(x, ..., oriented = FALSE, check = TRUE) { - new_s2_xptr( - s2_geography_from_wkb(x, oriented = oriented, check = check), - "s2_geography" - ) + s2_geog_from_wkb(x, oriented = oriented, check = check) } #' @rdname as_s2_geography #' @export as_s2_geography.blob <- function(x, ..., oriented = FALSE, check = TRUE) { - new_s2_xptr( - s2_geography_from_wkb(x, oriented = oriented, check = check), - "s2_geography" - ) + s2_geog_from_wkb(x, oriented = oriented, check = check) } #' @rdname as_s2_geography @@ -119,19 +113,16 @@ as_s2_geography.wk_wkt <- function(x, ..., oriented = FALSE, check = TRUE) { } } - new_s2_xptr( - s2_geography_from_wkt(x, oriented = oriented, check = check), - "s2_geography" + wk::wk_handle( + x, + s2_geography_writer(oriented = oriented, check = check) ) } #' @rdname as_s2_geography #' @export as_s2_geography.character <- function(x, ..., oriented = FALSE, check = TRUE) { - new_s2_xptr( - s2_geography_from_wkt(x, oriented = oriented, check = check), - "s2_geography" - ) + s2_geog_from_text(x, oriented = oriented, check = check) } #' @rdname as_s2_geography @@ -145,22 +136,20 @@ as_s2_geography.logical <- function(x, ...) { #' @rdname as_s2_geography #' @export as_wkb.s2_geography <- function(x, ...) { - wk::new_wk_wkb( - s2_geography_to_wkb(x, wk::wk_platform_endian()), - crs = wk::wk_crs_longlat("WGS84"), - geodesic = TRUE - ) + wkb <- wk::wk_handle(x, wk::wkb_writer()) + wk::wk_is_geodesic(wkb) <- TRUE + wk::wk_crs(wkb) <- wk::wk_crs_longlat() + wkb } #' @importFrom wk as_wkt #' @rdname as_s2_geography #' @export as_wkt.s2_geography <- function(x, ...) { - wk::new_wk_wkt( - s2_geography_to_wkt(x, precision = 16, trim = TRUE), - crs = wk::wk_crs_longlat(), - geodesic = TRUE - ) + wkt <- wk::wk_handle(x, wk::wkt_writer()) + wk::wk_is_geodesic(wkt) <- TRUE + wk::wk_crs(wkt) <- wk::wk_crs_longlat() + wkt } @@ -180,7 +169,7 @@ as_wkt.s2_geography <- function(x, ...) { #' @export format.s2_geography <- function(x, ..., max_coords = 5, precision = 9, trim = TRUE) { - paste0("<", s2_geography_format(x, max_coords, precision, trim), ">") + wk::wk_format(x, precision = precision, max_coords = max_coords, trim = trim) } # this is what gets called by the RStudio viewer, for which diff --git a/R/utils.R b/R/utils.R index 8d6d38e6..f12af485 100644 --- a/R/utils.R +++ b/R/utils.R @@ -68,8 +68,18 @@ stop_problems <- function(feature_id, problem, header) { expect_wkt_equal <- function(x, y, precision = 16) { testthat::expect_equal( - s2_geography_to_wkt(as_s2_geography(x), precision = precision, trim = TRUE), - s2_geography_to_wkt(as_s2_geography(y), precision = precision, trim = TRUE) + wk::wk_format( + as_s2_geography(x), + precision = precision, + trim = TRUE, + max_coords = .Machine$integer.max + ), + wk::wk_format( + as_s2_geography(y), + precision = precision, + trim = TRUE, + max_coords = .Machine$integer.max + ) ) } diff --git a/R/wk-utils.R b/R/wk-utils.R index 1372a361..8dda9bf6 100644 --- a/R/wk-utils.R +++ b/R/wk-utils.R @@ -65,6 +65,23 @@ s2_projection_filter <- function(handler, projection = s2_projection_plate_carre ) } +#' @importFrom wk wk_handle +#' @export +wk_handle.s2_geography <- function(geog, handler, ...) { + .Call(c_s2_handle_geography, geog, wk::as_wk_handler(handler)) +} + +s2_geography_writer <- function(oriented = FALSE, check = TRUE) { + wk::new_wk_handler( + .Call( + c_s2_geography_writer_new, + as.logical(oriented)[1], + as.logical(check)[1] + ), + "s2_geography_writer" + ) +} + #' @rdname s2_unprojection_filter #' @export s2_projection_plate_carree <- function() { diff --git a/data-raw/test-geom.R b/data-raw/test-geom.R new file mode 100644 index 00000000..b8b0c64f --- /dev/null +++ b/data-raw/test-geom.R @@ -0,0 +1,164 @@ + +nc_wkt <- wk::as_wkt(sf::read_sf(system.file("shape/nc.shp", package = "sf"))) +wk::wk_crs(nc_wkt) <- wk::wk_crs_proj_definition(wk::wk_crs(nc_wkt), verbose = TRUE) + +geoarrow_example_wkt <- list( + nc = nc_wkt, + point = c("POINT (30 10)", "POINT EMPTY", NA), + linestring = c( + "LINESTRING (30 10, 10 30, 40 40)", "LINESTRING EMPTY", + NA + ), + polygon = c( + "POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))", + "POLYGON ((35 10, 45 45, 15 40, 10 20, 35 10), (20 30, 35 35, 30 20, 20 30))", + "POLYGON EMPTY", NA + ), + multipoint = c( + "MULTIPOINT ((30 10))", + "MULTIPOINT ((10 40), (40 30), (20 20), (30 10))", "MULTIPOINT ((10 40), (40 30), (20 20), (30 10))", + "MULTIPOINT EMPTY", NA + ), + multilinestring = c( + "MULTILINESTRING ((30 10, 10 30, 40 40))", + "MULTILINESTRING ((10 10, 20 20, 10 40), (40 40, 30 30, 40 20, 30 10))", + "MULTILINESTRING EMPTY", NA + ), + multipolygon = c( + "MULTIPOLYGON (((30 10, 40 40, 20 40, 10 20, 30 10)))", + "MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)), ((15 5, 40 10, 10 20, 5 10, 15 5)))", + "MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)), ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))", + "MULTIPOLYGON EMPTY", NA + ), + geometrycollection = c( + "GEOMETRYCOLLECTION (POINT (30 10))", + "GEOMETRYCOLLECTION (LINESTRING (30 10, 10 30, 40 40))", "GEOMETRYCOLLECTION (POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10)))", + "GEOMETRYCOLLECTION (POINT (30 10), LINESTRING (30 10, 10 30, 40 40), POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10)))", + "GEOMETRYCOLLECTION (GEOMETRYCOLLECTION (POINT (30 10)))", "GEOMETRYCOLLECTION (GEOMETRYCOLLECTION (LINESTRING (30 10, 10 30, 40 40)))", + "GEOMETRYCOLLECTION (GEOMETRYCOLLECTION (POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))))", + "GEOMETRYCOLLECTION EMPTY", NA + ), + point_z = c( + "POINT Z (30 10 40)", + "POINT Z EMPTY", NA + ), + linestring_z = c( + "LINESTRING Z (30 10 40, 10 30 40, 40 40 80)", + "LINESTRING EMPTY", NA + ), + polygon_z = c( + "POLYGON Z ((30 10 40, 40 40 80, 20 40 60, 10 20 30, 30 10 40))", + "POLYGON Z ((35 10 45, 45 45 90, 15 40 55, 10 20 30, 35 10 45), (20 30 50, 35 35 70, 30 20 50, 20 30 50))", + "POLYGON EMPTY", NA + ), + multipoint_z = c( + "MULTIPOINT Z ((30 10 40))", + "MULTIPOINT Z ((10 40 50), (40 30 70), (20 20 40), (30 10 40))", + "MULTIPOINT Z ((10 40 50), (40 30 70), (20 20 40), (30 10 40))", + "MULTIPOINT EMPTY", NA + ), + multilinestring_z = c( + "MULTILINESTRING Z ((30 10 40, 10 30 40, 40 40 80))", + "MULTILINESTRING Z ((10 10 20, 20 20 40, 10 40 50), (40 40 80, 30 30 60, 40 20 60, 30 10 40))", + "MULTILINESTRING EMPTY", NA + ), + multipolygon_z = c( + "MULTIPOLYGON Z (((30 10 40, 40 40 80, 20 40 60, 10 20 30, 30 10 40)))", + "MULTIPOLYGON Z (((30 20 50, 45 40 85, 10 40 50, 30 20 50)), ((15 5 20, 40 10 50, 10 20 30, 5 10 15, 15 5 20)))", + "MULTIPOLYGON Z (((40 40 80, 20 45 65, 45 30 75, 40 40 80)), ((20 35 55, 10 30 40, 10 10 20, 30 5 35, 45 20 65, 20 35 55), (30 20 50, 20 15 35, 20 25 45, 30 20 50)))", + "MULTIPOLYGON EMPTY", NA + ), + geometrycollection_z = c( + "GEOMETRYCOLLECTION Z (POINT Z (30 10 40))", + "GEOMETRYCOLLECTION Z (LINESTRING Z (30 10 40, 10 30 40, 40 40 80))", + "GEOMETRYCOLLECTION Z (POLYGON Z ((30 10 40, 40 40 80, 20 40 60, 10 20 30, 30 10 40)))", + "GEOMETRYCOLLECTION Z (POINT Z (30 10 40), LINESTRING Z (30 10 40, 10 30 40, 40 40 80), POLYGON Z ((30 10 40, 40 40 80, 20 40 60, 10 20 30, 30 10 40)))", + "GEOMETRYCOLLECTION Z (GEOMETRYCOLLECTION Z (POINT Z (30 10 40)))", + "GEOMETRYCOLLECTION Z (GEOMETRYCOLLECTION Z (LINESTRING Z (30 10 40, 10 30 40, 40 40 80)))", + "GEOMETRYCOLLECTION Z (GEOMETRYCOLLECTION Z (POLYGON Z ((30 10 40, 40 40 80, 20 40 60, 10 20 30, 30 10 40))))", + "GEOMETRYCOLLECTION EMPTY", NA + ), + point_m = c( + "POINT M (30 10 300)", + "POINT M EMPTY", NA + ), linestring_m = c( + "LINESTRING M (30 10 300, 10 30 300, 40 40 1600)", + "LINESTRING EMPTY", NA + ), + polygon_m = c( + "POLYGON M ((30 10 300, 40 40 1600, 20 40 800, 10 20 200, 30 10 300))", + "POLYGON M ((35 10 350, 45 45 2025, 15 40 600, 10 20 200, 35 10 350), (20 30 600, 35 35 1225, 30 20 600, 20 30 600))", + "POLYGON EMPTY", NA + ), + multipoint_m = c( + "MULTIPOINT M ((30 10 300))", + "MULTIPOINT M ((10 40 400), (40 30 1200), (20 20 400), (30 10 300))", + "MULTIPOINT M ((10 40 400), (40 30 1200), (20 20 400), (30 10 300))", + "MULTIPOINT EMPTY", NA + ), + multilinestring_m = c( + "MULTILINESTRING M ((30 10 300, 10 30 300, 40 40 1600))", + "MULTILINESTRING M ((10 10 100, 20 20 400, 10 40 400), (40 40 1600, 30 30 900, 40 20 800, 30 10 300))", + "MULTILINESTRING EMPTY", NA + ), + multipolygon_m = c( + "MULTIPOLYGON M (((30 10 300, 40 40 1600, 20 40 800, 10 20 200, 30 10 300)))", + "MULTIPOLYGON M (((30 20 600, 45 40 1800, 10 40 400, 30 20 600)), ((15 5 75, 40 10 400, 10 20 200, 5 10 50, 15 5 75)))", + "MULTIPOLYGON M (((40 40 1600, 20 45 900, 45 30 1350, 40 40 1600)), ((20 35 700, 10 30 300, 10 10 100, 30 5 150, 45 20 900, 20 35 700), (30 20 600, 20 15 300, 20 25 500, 30 20 600)))", + "MULTIPOLYGON EMPTY", NA + ), + geometrycollection_m = c( + "GEOMETRYCOLLECTION M (POINT M (30 10 300))", + "GEOMETRYCOLLECTION M (LINESTRING M (30 10 300, 10 30 300, 40 40 1600))", + "GEOMETRYCOLLECTION M (POLYGON M ((30 10 300, 40 40 1600, 20 40 800, 10 20 200, 30 10 300)))", + "GEOMETRYCOLLECTION M (POINT M (30 10 300), LINESTRING M (30 10 300, 10 30 300, 40 40 1600), POLYGON M ((30 10 300, 40 40 1600, 20 40 800, 10 20 200, 30 10 300)))", + "GEOMETRYCOLLECTION M (GEOMETRYCOLLECTION M (POINT M (30 10 300)))", + "GEOMETRYCOLLECTION M (GEOMETRYCOLLECTION M (LINESTRING M (30 10 300, 10 30 300, 40 40 1600)))", + "GEOMETRYCOLLECTION M (GEOMETRYCOLLECTION M (POLYGON M ((30 10 300, 40 40 1600, 20 40 800, 10 20 200, 30 10 300))))", + "GEOMETRYCOLLECTION EMPTY", NA + ), + point_zm = c( + "POINT ZM (30 10 40 300)", + "POINT ZM EMPTY", NA + ), + linestring_zm = c( + "LINESTRING ZM (30 10 40 300, 10 30 40 300, 40 40 80 1600)", + "LINESTRING EMPTY", NA + ), + polygon_zm = c( + "POLYGON ZM ((30 10 40 300, 40 40 80 1600, 20 40 60 800, 10 20 30 200, 30 10 40 300))", + "POLYGON ZM ((35 10 45 350, 45 45 90 2025, 15 40 55 600, 10 20 30 200, 35 10 45 350), (20 30 50 600, 35 35 70 1225, 30 20 50 600, 20 30 50 600))", + "POLYGON EMPTY", NA + ), + multipoint_zm = c( + "MULTIPOINT ZM ((30 10 40 300))", + "MULTIPOINT ZM ((10 40 50 400), (40 30 70 1200), (20 20 40 400), (30 10 40 300))", + "MULTIPOINT ZM ((10 40 50 400), (40 30 70 1200), (20 20 40 400), (30 10 40 300))", + "MULTIPOINT EMPTY", NA + ), + multilinestring_zm = c( + "MULTILINESTRING ZM ((30 10 40 300, 10 30 40 300, 40 40 80 1600))", + "MULTILINESTRING ZM ((10 10 20 100, 20 20 40 400, 10 40 50 400), (40 40 80 1600, 30 30 60 900, 40 20 60 800, 30 10 40 300))", + "MULTILINESTRING EMPTY", NA + ), + multipolygon_zm = c( + "MULTIPOLYGON ZM (((30 10 40 300, 40 40 80 1600, 20 40 60 800, 10 20 30 200, 30 10 40 300)))", + "MULTIPOLYGON ZM (((30 20 50 600, 45 40 85 1800, 10 40 50 400, 30 20 50 600)), ((15 5 20 75, 40 10 50 400, 10 20 30 200, 5 10 15 50, 15 5 20 75)))", + "MULTIPOLYGON ZM (((40 40 80 1600, 20 45 65 900, 45 30 75 1350, 40 40 80 1600)), ((20 35 55 700, 10 30 40 300, 10 10 20 100, 30 5 35 150, 45 20 65 900, 20 35 55 700), (30 20 50 600, 20 15 35 300, 20 25 45 500, 30 20 50 600)))", + "MULTIPOLYGON EMPTY", NA + ), + geometrycollection_zm = c( + "GEOMETRYCOLLECTION ZM (POINT ZM (30 10 40 300))", + "GEOMETRYCOLLECTION ZM (LINESTRING ZM (30 10 40 300, 10 30 40 300, 40 40 80 1600))", + "GEOMETRYCOLLECTION ZM (POLYGON ZM ((30 10 40 300, 40 40 80 1600, 20 40 60 800, 10 20 30 200, 30 10 40 300)))", + "GEOMETRYCOLLECTION ZM (POINT ZM (30 10 40 300), LINESTRING ZM (30 10 40 300, 10 30 40 300, 40 40 80 1600), POLYGON ZM ((30 10 40 300, 40 40 80 1600, 20 40 60 800, 10 20 30 200, 30 10 40 300)))", + "GEOMETRYCOLLECTION ZM (GEOMETRYCOLLECTION ZM (POINT ZM (30 10 40 300)))", + "GEOMETRYCOLLECTION ZM (GEOMETRYCOLLECTION ZM (LINESTRING ZM (30 10 40 300, 10 30 40 300, 40 40 80 1600)))", + "GEOMETRYCOLLECTION ZM (GEOMETRYCOLLECTION ZM (POLYGON ZM ((30 10 40 300, 40 40 80 1600, 20 40 60 800, 10 20 30 200, 30 10 40 300))))", + "GEOMETRYCOLLECTION EMPTY", NA + ) +) + +s2_data_example_wkt <- lapply(geoarrow_example_wkt, wk::wkt, geodesic = TRUE) + +usethis::use_data(s2_data_example_wkt, overwrite = TRUE) diff --git a/data/s2_data_example_wkt.rda b/data/s2_data_example_wkt.rda new file mode 100644 index 00000000..23820bbf Binary files /dev/null and b/data/s2_data_example_wkt.rda differ diff --git a/man/s2_data_example_wkt.Rd b/man/s2_data_example_wkt.Rd new file mode 100644 index 00000000..1dea278c --- /dev/null +++ b/man/s2_data_example_wkt.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{s2_data_example_wkt} +\alias{s2_data_example_wkt} +\title{Example Geometries} +\format{ +An object of class \code{list} of length 29. +} +\usage{ +s2_data_example_wkt +} +\description{ +These geometries are toy examples useful for testing various coordinate +shuffling operations in the s2 package. +} +\keyword{datasets} diff --git a/src/Makevars.in b/src/Makevars.in index 943c0de8..672304cb 100644 --- a/src/Makevars.in +++ b/src/Makevars.in @@ -103,7 +103,15 @@ OBJECTS = $(ABSL_LIBS) \ s2-xptr.o \ wk-impl.o \ wk-c-utils.o \ - s2-c-api.o \ + tessellator.o \ + s2-geography/accessors.o \ + s2-geography/accessors-geog.o \ + s2-geography/linear-referencing.o \ + s2-geography/distance.o \ + s2-geography/build.o \ + s2-geography/coverings.o \ + s2-geography/geography.o \ + s2-geography/predicates.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 e87c2d08..15a8af20 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -88,7 +88,14 @@ ABSL_LIBS = absl/base/internal/cycleclock.o \ absl/types/bad_variant_access.o S2LIBS = $(ABSL_LIBS) \ - s2-c-api.o \ + s2-geography/linear-referencing.o \ + s2-geography/distance.o \ + s2-geography/accessors.o \ + s2-geography/accessors-geog.o \ + s2-geography/build.o \ + s2-geography/coverings.o \ + s2-geography/geography.o \ + s2-geography/predicates.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 c018de11..cd87e86d 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -737,73 +737,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// cpp_s2_geog_point -List cpp_s2_geog_point(NumericVector x, NumericVector y); -RcppExport SEXP _s2_cpp_s2_geog_point(SEXP xSEXP, SEXP ySEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - rcpp_result_gen = Rcpp::wrap(cpp_s2_geog_point(x, y)); - return rcpp_result_gen; -END_RCPP -} -// cpp_s2_make_line -List cpp_s2_make_line(NumericVector x, NumericVector y, IntegerVector featureId); -RcppExport SEXP _s2_cpp_s2_make_line(SEXP xSEXP, SEXP ySEXP, SEXP featureIdSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - Rcpp::traits::input_parameter< IntegerVector >::type featureId(featureIdSEXP); - rcpp_result_gen = Rcpp::wrap(cpp_s2_make_line(x, y, featureId)); - return rcpp_result_gen; -END_RCPP -} -// cpp_s2_make_polygon -List cpp_s2_make_polygon(NumericVector x, NumericVector y, IntegerVector featureId, IntegerVector ringId, bool oriented, bool check); -RcppExport SEXP _s2_cpp_s2_make_polygon(SEXP xSEXP, SEXP ySEXP, SEXP featureIdSEXP, SEXP ringIdSEXP, SEXP orientedSEXP, SEXP checkSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); - Rcpp::traits::input_parameter< NumericVector >::type y(ySEXP); - Rcpp::traits::input_parameter< IntegerVector >::type featureId(featureIdSEXP); - Rcpp::traits::input_parameter< IntegerVector >::type ringId(ringIdSEXP); - Rcpp::traits::input_parameter< bool >::type oriented(orientedSEXP); - Rcpp::traits::input_parameter< bool >::type check(checkSEXP); - rcpp_result_gen = Rcpp::wrap(cpp_s2_make_polygon(x, y, featureId, ringId, oriented, check)); - return rcpp_result_gen; -END_RCPP -} -// s2_geography_from_wkb -List s2_geography_from_wkb(List wkb, bool oriented, bool check); -RcppExport SEXP _s2_s2_geography_from_wkb(SEXP wkbSEXP, SEXP orientedSEXP, SEXP checkSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< List >::type wkb(wkbSEXP); - Rcpp::traits::input_parameter< bool >::type oriented(orientedSEXP); - Rcpp::traits::input_parameter< bool >::type check(checkSEXP); - rcpp_result_gen = Rcpp::wrap(s2_geography_from_wkb(wkb, oriented, check)); - return rcpp_result_gen; -END_RCPP -} -// s2_geography_from_wkt -List s2_geography_from_wkt(CharacterVector wkt, bool oriented, bool check); -RcppExport SEXP _s2_s2_geography_from_wkt(SEXP wktSEXP, SEXP orientedSEXP, SEXP checkSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< CharacterVector >::type wkt(wktSEXP); - Rcpp::traits::input_parameter< bool >::type oriented(orientedSEXP); - Rcpp::traits::input_parameter< bool >::type check(checkSEXP); - rcpp_result_gen = Rcpp::wrap(s2_geography_from_wkt(wkt, oriented, check)); - return rcpp_result_gen; -END_RCPP -} // s2_geography_full List s2_geography_full(LogicalVector x); RcppExport SEXP _s2_s2_geography_full(SEXP xSEXP) { @@ -815,45 +748,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// s2_geography_to_wkt -CharacterVector s2_geography_to_wkt(List s2_geography, int precision, bool trim); -RcppExport SEXP _s2_s2_geography_to_wkt(SEXP s2_geographySEXP, SEXP precisionSEXP, SEXP trimSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< List >::type s2_geography(s2_geographySEXP); - Rcpp::traits::input_parameter< int >::type precision(precisionSEXP); - Rcpp::traits::input_parameter< bool >::type trim(trimSEXP); - rcpp_result_gen = Rcpp::wrap(s2_geography_to_wkt(s2_geography, precision, trim)); - return rcpp_result_gen; -END_RCPP -} -// s2_geography_to_wkb -List s2_geography_to_wkb(List s2_geography, int endian); -RcppExport SEXP _s2_s2_geography_to_wkb(SEXP s2_geographySEXP, SEXP endianSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< List >::type s2_geography(s2_geographySEXP); - Rcpp::traits::input_parameter< int >::type endian(endianSEXP); - rcpp_result_gen = Rcpp::wrap(s2_geography_to_wkb(s2_geography, endian)); - return rcpp_result_gen; -END_RCPP -} -// s2_geography_format -CharacterVector s2_geography_format(List s2_geography, int maxCoords, int precision, bool trim); -RcppExport SEXP _s2_s2_geography_format(SEXP s2_geographySEXP, SEXP maxCoordsSEXP, SEXP precisionSEXP, SEXP trimSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< List >::type s2_geography(s2_geographySEXP); - Rcpp::traits::input_parameter< int >::type maxCoords(maxCoordsSEXP); - Rcpp::traits::input_parameter< int >::type precision(precisionSEXP); - Rcpp::traits::input_parameter< bool >::type trim(trimSEXP); - rcpp_result_gen = Rcpp::wrap(s2_geography_format(s2_geography, maxCoords, precision, trim)); - return rcpp_result_gen; -END_RCPP -} // s2_lnglat_from_numeric List s2_lnglat_from_numeric(NumericVector lng, NumericVector lat); RcppExport SEXP _s2_s2_lnglat_from_numeric(SEXP lngSEXP, SEXP latSEXP) { @@ -1482,6 +1376,8 @@ END_RCPP } RcppExport SEXP c_s2_coord_filter_new(SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP c_s2_geography_writer_new(SEXP, SEXP); +RcppExport SEXP c_s2_handle_geography(SEXP, SEXP); RcppExport SEXP c_s2_projection_mercator(); RcppExport SEXP c_s2_projection_plate_carree(); @@ -1549,15 +1445,7 @@ static const R_CallMethodDef CallEntries[] = { {"_s2_cpp_s2_cell_max_distance", (DL_FUNC) &_s2_cpp_s2_cell_max_distance, 2}, {"_s2_cpp_s2_cell_common_ancestor_level", (DL_FUNC) &_s2_cpp_s2_cell_common_ancestor_level, 2}, {"_s2_cpp_s2_cell_common_ancestor_level_agg", (DL_FUNC) &_s2_cpp_s2_cell_common_ancestor_level_agg, 1}, - {"_s2_cpp_s2_geog_point", (DL_FUNC) &_s2_cpp_s2_geog_point, 2}, - {"_s2_cpp_s2_make_line", (DL_FUNC) &_s2_cpp_s2_make_line, 3}, - {"_s2_cpp_s2_make_polygon", (DL_FUNC) &_s2_cpp_s2_make_polygon, 6}, - {"_s2_s2_geography_from_wkb", (DL_FUNC) &_s2_s2_geography_from_wkb, 3}, - {"_s2_s2_geography_from_wkt", (DL_FUNC) &_s2_s2_geography_from_wkt, 3}, {"_s2_s2_geography_full", (DL_FUNC) &_s2_s2_geography_full, 1}, - {"_s2_s2_geography_to_wkt", (DL_FUNC) &_s2_s2_geography_to_wkt, 3}, - {"_s2_s2_geography_to_wkb", (DL_FUNC) &_s2_s2_geography_to_wkb, 2}, - {"_s2_s2_geography_format", (DL_FUNC) &_s2_s2_geography_format, 4}, {"_s2_s2_lnglat_from_numeric", (DL_FUNC) &_s2_s2_lnglat_from_numeric, 2}, {"_s2_s2_lnglat_from_s2_point", (DL_FUNC) &_s2_s2_lnglat_from_s2_point, 1}, {"_s2_data_frame_from_s2_lnglat", (DL_FUNC) &_s2_data_frame_from_s2_lnglat, 1}, @@ -1609,6 +1497,8 @@ 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_geography_writer_new", (DL_FUNC) &c_s2_geography_writer_new, 2}, + {"c_s2_handle_geography", (DL_FUNC) &c_s2_handle_geography, 2}, {"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/geography-collection.h b/src/geography-collection.h deleted file mode 100644 index 3c2b6cbd..00000000 --- a/src/geography-collection.h +++ /dev/null @@ -1,245 +0,0 @@ - -#ifndef GEOGRAPHY_COLLECTION_H -#define GEOGRAPHY_COLLECTION_H - -#include -#include "geography.h" - -// This class handles collections of other Geography -// objects. -class GeographyCollection: public Geography { -public: - GeographyCollection(): features(0) {} - GeographyCollection(std::vector> features): - features(std::move(features)) {} - - Geography::Type GeographyType() { - return Geography::Type::GEOGRAPHY_COLLECTION; - } - - bool FindValidationError(S2Error* error) { - bool result; - error->Clear(); - for (size_t i = 0; i < this->features.size(); i++) { - result = this->features[i]->FindValidationError(error); - if (result) { - return result; - } - } - - return false; - } - - bool IsCollection() { - return this->features.size() > 0; - } - - int Dimension() { - int dimension = -1; - for (size_t i = 0; i < this->features.size(); i++) { - dimension = std::max(this->features[i]->Dimension(), dimension); - } - - return dimension; - } - - int NumPoints() { - int numPoints = 0; - for (size_t i = 0; i < this->features.size(); i++) { - numPoints += this->features[i]->NumPoints(); - } - return numPoints; - } - - bool IsEmpty() { - // could also loop and test all(!IsEmpty()), but - // that is inconsistent with what gets printed - return this->features.size() == 0; - } - - double Area() { - double area = 0; - for (size_t i = 0; i < this->features.size(); i++) { - area += this->features[i]->Area(); - } - return area; - } - - double Length() { - double length = 0; - for (size_t i = 0; i < this->features.size(); i++) { - length += this->features[i]->Length(); - } - return length; - } - - double Perimeter() { - double perimeter = 0; - for (size_t i = 0; i < this->features.size(); i++) { - perimeter += this->features[i]->Perimeter(); - } - return perimeter; - } - - double X() { - Rcpp::stop("Can't compute X value of a non-point geography"); - } - - double Y() { - Rcpp::stop("Can't compute Y value of a non-point geography"); - } - - S2Point Centroid() { - S2Point cumCentroid(0, 0, 0); - for (size_t i = 0; i < this->features.size(); i++) { - S2Point centroid = this->features[i]->Centroid(); - if (centroid.Norm2() > 0) { - cumCentroid += centroid.Normalize(); - } - } - - return cumCentroid; - } - - std::unique_ptr Boundary() { - std::vector> featureBoundaries(this->features.size()); - for (size_t i = 0; i < this->features.size(); i++) { - featureBoundaries[i] = this->features[i]->Boundary(); - } - - return absl::make_unique(std::move(featureBoundaries)); - } - - std::vector BuildShapeIndex(MutableS2ShapeIndex* index) { - std::vector shapeIds; - for (size_t i = 0; i < this->features.size(); i++) { - std::vector newShapeIds = this->features[i]->BuildShapeIndex(index); - for (size_t j = 0; j < newShapeIds.size(); j++) { - shapeIds.push_back(newShapeIds[j]); - } - } - return shapeIds; - } - - const std::vector >* CollectionFeatures() { - return &(this->features); - } - - void Export(WKGeometryHandler* handler, uint32_t partId) { - WKGeometryMeta meta(WKGeometryType::GeometryCollection, false, false, false); - meta.hasSize = true; - meta.size = this->features.size(); - - handler->nextGeometryStart(meta, partId); - for (size_t i = 0; i < this->features.size(); i++) { - this->features[i]->Export(handler, i); - } - handler->nextGeometryEnd(meta, partId); - } - - class Builder: public GeographyBuilder { - public: - - Builder(bool oriented, bool check): - metaPtr(nullptr), builderPtr(nullptr), builderMetaPtr(nullptr), - oriented(oriented), check(check) {} - - virtual void nextGeometryStart(const WKGeometryMeta& meta, uint32_t partId) { - // if this is the first call, store the meta reference associated with this geometry - if (this->metaPtr == nullptr) { - this->metaPtr = (WKGeometryMeta*) &meta; - return; - } - - if (!this->builderPtr) { - // store a reference to the meta associated with this - // builder so that we know when the corresponding nextGeometryEnd() - // is called - this->builderMetaPtr = (WKGeometryMeta*) &meta; - - switch (meta.geometryType) { - case WKGeometryType::Point: - case WKGeometryType::MultiPoint: - this->builderPtr = absl::make_unique(); - break; - case WKGeometryType::LineString: - case WKGeometryType::MultiLineString: - this->builderPtr = absl::make_unique(); - break; - case WKGeometryType::Polygon: - case WKGeometryType::MultiPolygon: - this->builderPtr = absl::make_unique( - this->oriented, - this->check - ); - break; - case WKGeometryType::GeometryCollection: - this->builderPtr = absl::make_unique( - this->oriented, - this->check - ); - break; - default: - std::stringstream err; - err << "Unknown geometry type in geography builder: " << meta.geometryType; - Rcpp::stop(err.str()); - } - } - - this->builder()->nextGeometryStart(meta, partId); - } - - virtual void nextLinearRingStart(const WKGeometryMeta& meta, uint32_t size, uint32_t ringId) { - this->builder()->nextLinearRingStart(meta, size, ringId); - } - - virtual void nextCoordinate(const WKGeometryMeta& meta, const WKCoord& coord, uint32_t coordId) { - this->builder()->nextCoordinate(meta, coord, coordId); - } - - virtual void nextLinearRingEnd(const WKGeometryMeta& meta, uint32_t size, uint32_t ringId) { - this->builder()->nextLinearRingEnd(meta, size, ringId); - } - - virtual void nextGeometryEnd(const WKGeometryMeta& meta, uint32_t partId) { - // the end of this GEOMETRYCOLLECTION - if (&meta == this->metaPtr) { - return; - } - - this->builder()->nextGeometryEnd(meta, partId); - - if (&meta == this->builderMetaPtr) { - std::unique_ptr feature = this->builder()->build(); - features.push_back(std::move(feature)); - this->builderPtr = std::unique_ptr(nullptr); - this->builderMetaPtr = nullptr; - } - } - - std::unique_ptr build() { - return absl::make_unique(std::move(this->features)); - } - - private: - std::vector> features; - WKGeometryMeta* metaPtr; - std::unique_ptr builderPtr; - WKGeometryMeta* builderMetaPtr; - bool oriented; - bool check; - - GeographyBuilder* builder() { - if (this->builderPtr) { - return this->builderPtr.get(); - } else { - Rcpp::stop("Invalid nesting in geometrycollection (can't find nested builder)"); - } - } - }; - -private: - std::vector> features; -}; - -#endif diff --git a/src/geography.h b/src/geography.h index 20f47149..a20f7394 100644 --- a/src/geography.h +++ b/src/geography.h @@ -2,119 +2,94 @@ #ifndef GEOGRAPHY_H #define GEOGRAPHY_H -#include -#include "s2/s2latlng.h" -#include "s2/s2polyline.h" -#include "s2/s2polygon.h" -#include "s2/s2shape_index.h" -#include "s2/s2shape_index_region.h" -#include "s2/mutable_s2shape_index.h" -#include "s2/s2point_vector_shape.h" -#include "s2/s2cap.h" -#include "wk/geometry-handler.hpp" -#include +#define R_NO_REMAP +#include +#include + +#include "s2-geography/s2-geography.hpp" class Geography { public: + Geography(std::unique_ptr geog): + geog_(std::move(geog)), index_(nullptr) {} - enum class Type { - GEOGRAPHY_EMPTY, - GEOGRAPHY_POINT, - GEOGRAPHY_POLYLINE, - GEOGRAPHY_POLYGON, - GEOGRAPHY_COLLECTION - }; + const s2geography::S2Geography& Geog() const { + return *geog_; + } - Geography(): hasIndex(false) {} + const s2geography::ShapeIndexGeography& Index() { + if (!index_) { + this->index_ = absl::make_unique(*geog_); + } - // accessors need to be methods, since their calculation - // depends on the geometry type - virtual Type GeographyType() { - return Type::GEOGRAPHY_EMPTY; + return *index_; } - virtual bool FindValidationError(S2Error* error) = 0; - - // returns true for a multi- - // or geometrycollection type - virtual bool IsCollection() = 0; - // Returns 0 for point, 1 for line, 2 for polygon - virtual int Dimension() = 0; - // Returns the number of points in the input - virtual int NumPoints() = 0; - virtual bool IsEmpty() = 0; - virtual double Area() = 0; - virtual double Length() = 0; - virtual double Perimeter() = 0; - virtual double X() = 0; - virtual double Y() = 0; - virtual S2Point Centroid() = 0; - virtual std::unique_ptr Boundary() = 0; - - // every type will build the index differently based on - // the underlying data, and this can (should?) be done - // lazily. Returns a vector of shape IDs so the caller - // can keep track of which shape came from which feature. - virtual std::vector BuildShapeIndex(MutableS2ShapeIndex* index) = 0; - - // the factory handler is responsible for building these objects - // but exporting can be done here - virtual void Export(WKGeometryHandler* handler, uint32_t partId) = 0; - - virtual ~Geography() {} - - // Most calculations will use the ShapeIndex, but sometimes access to the - // underlying point, line, or polygon is useful to keep this class from - // becoming bloated with the entire s2 API. - virtual const std::vector* Point() { - return nullptr; + static SEXP MakeXPtr(std::unique_ptr geog) { + SEXP xptr = PROTECT(R_MakeExternalPtr(new Geography(std::move(geog)), R_NilValue, R_NilValue)); + R_RegisterCFinalizer(xptr, &finalize_xptr); + UNPROTECT(1); + return xptr; } - virtual const std::vector>* Polyline() { - return nullptr; + static SEXP MakeXPtr(std::unique_ptr geog) { + std::unique_ptr geog_owning = std::move(geog); + SEXP xptr = PROTECT(R_MakeExternalPtr(geog_owning.release(), R_NilValue, R_NilValue)); + R_RegisterCFinalizer(xptr, &finalize_xptr); + UNPROTECT(1); + return xptr; } - virtual const S2Polygon* Polygon() { - return nullptr; + static std::unique_ptr MakePoint() { + return absl::make_unique(absl::make_unique()); } - virtual const std::vector>* CollectionFeatures() { - return nullptr; + static std::unique_ptr MakePoint(S2Point point) { + return absl::make_unique(absl::make_unique(point)); } - // other calculations use ShapeIndex - virtual S2ShapeIndex* ShapeIndex() { - if (!this->hasIndex) { - this->BuildShapeIndex(&this->shape_index_); - this->hasIndex = true; - } + static std::unique_ptr MakePoint(std::vector points) { + return absl::make_unique(absl::make_unique(std::move(points))); + } - return &this->shape_index_; + static std::unique_ptr MakePolyline() { + return absl::make_unique(absl::make_unique()); } - virtual S2ShapeIndexRegion ShapeIndexRegion() { - S2ShapeIndex *ix = this->ShapeIndex(); - return MakeS2ShapeIndexRegion(ix); + static std::unique_ptr MakePolyline(std::unique_ptr polyline) { + return absl::make_unique(absl::make_unique(std::move(polyline))); } - virtual S2Cap GetCapBound() { - return this->ShapeIndexRegion().GetCapBound(); + static std::unique_ptr MakePolyline(std::vector> polylines) { + return absl::make_unique(absl::make_unique(std::move(polylines))); } - virtual S2LatLngRect GetRectBound() { - return this->ShapeIndexRegion().GetRectBound(); + static std::unique_ptr MakePolygon() { + return absl::make_unique(absl::make_unique()); } -protected: - MutableS2ShapeIndex shape_index_; - bool hasIndex; -}; + static std::unique_ptr MakePolygon(std::unique_ptr polygon) { + return absl::make_unique(absl::make_unique(std::move(polygon))); + } + static std::unique_ptr MakeCollection() { + return absl::make_unique(absl::make_unique()); + } -class GeographyBuilder: public WKGeometryHandler { -public: - virtual std::unique_ptr build() = 0; - virtual ~GeographyBuilder() {} + static std::unique_ptr MakeCollection(std::vector> features) { + return absl::make_unique(absl::make_unique(std::move(features))); + } + +private: + std::unique_ptr geog_; + std::unique_ptr index_; + + static void finalize_xptr(SEXP xptr) { + Geography* geog = reinterpret_cast(R_ExternalPtrAddr(xptr)); + if (geog != nullptr) { + delete geog; + } + } }; #endif diff --git a/src/point-geography.h b/src/point-geography.h deleted file mode 100644 index c4f16f0a..00000000 --- a/src/point-geography.h +++ /dev/null @@ -1,181 +0,0 @@ - -#ifndef POINT_GEOGRAPHY_H -#define POINT_GEOGRAPHY_H - -#include - -#include "s2/s2latlng_rect.h" - -#include "geography.h" - -// This class handles both points and multipoints, as this is how -// points are generally returned/required in S2 (vector of S2Point) -// This is similar to an S2PointVectorLayer -class PointGeography: public Geography { -public: - PointGeography(): points(0) {} - PointGeography(S2Point point): points(1) { - this->points[0] = point; - } - PointGeography(std::vector points): points(points) {} - - Geography::Type GeographyType() { - return Geography::Type::GEOGRAPHY_POINT; - } - - bool FindValidationError(S2Error* error) { - return false; - } - - const std::vector* Point() { - return &(this->points); - } - - bool IsCollection() { - return this->points.size() > 1; - } - - int Dimension() { - return 0; - } - - int NumPoints() { - return this->points.size(); - } - - bool IsEmpty() { - return this->points.size() == 0; - } - - double Area() { - return 0; - } - - double Length() { - return 0; - } - - double Perimeter() { - return 0; - } - - double X() { - if (this->points.size() != 1) { - return NA_REAL; - } else { - S2LatLng latLng(this->points[0]); - return latLng.lng().degrees(); - } - } - - double Y() { - if (this->points.size() != 1) { - return NA_REAL; - } else { - S2LatLng latLng(this->points[0]); - return latLng.lat().degrees(); - } - } - - S2Point Centroid() { - S2Point output(0, 0, 0); - for (size_t i = 0; i < this->points.size(); i++) { - output += this->points[i]; - } - - return output; - } - - S2LatLngRect GetRectBound() { - S2LatLngRect rect; - for (size_t i = 0; i < this->points.size(); i++) { - rect.AddPoint(this->points[i]); // depends on order - } - return rect; - } - - std::unique_ptr Boundary() { - return absl::make_unique(); - } - - std::vector BuildShapeIndex(MutableS2ShapeIndex* index) { - std::vector shapeIds(1); - std::vector pointsCopy(this->points); - - shapeIds[0] = index->Add(std::unique_ptr( - new S2PointVectorShape(std::move(pointsCopy))) - ); - return shapeIds; - } - - void Export(WKGeometryHandler* handler, uint32_t partId) { - S2LatLng point; - - if (this->points.size() > 1) { - // export multipoint - WKGeometryMeta meta(WKGeometryType::MultiPoint, false, false, false); - meta.hasSize = true; - meta.size = this->points.size(); - - WKGeometryMeta childMeta(WKGeometryType::Point, false, false, false); - childMeta.hasSize = true; - childMeta.size = 1; - - handler->nextGeometryStart(meta, partId); - - for (size_t i = 0; i < this->points.size(); i++) { - point = S2LatLng(this->points[i]); - - handler->nextGeometryStart(childMeta, i); - handler->nextCoordinate(meta, WKCoord::xy(point.lng().degrees(), point.lat().degrees()), 0); - handler->nextGeometryEnd(childMeta, i); - } - - handler->nextGeometryEnd(meta, partId); - - } else if (this->points.size() > 0) { - // export point - WKGeometryMeta meta(WKGeometryType::Point, false, false, false); - meta.hasSize = true; - meta.size = this->points.size(); - - handler->nextGeometryStart(meta, partId); - - point = S2LatLng(this->points[0]); - handler->nextCoordinate(meta, WKCoord::xy(point.lng().degrees(), point.lat().degrees()), 0); - - handler->nextGeometryEnd(meta, partId); - } else { - // export empty point - // export point - WKGeometryMeta meta(WKGeometryType::Point, false, false, false); - meta.hasSize = true; - meta.size = 0; - handler->nextGeometryStart(meta, partId); - handler->nextGeometryEnd(meta, partId); - } - } - - class Builder: public GeographyBuilder { - public: - void nextCoordinate(const WKGeometryMeta& meta, const WKCoord& coord, uint32_t coordId) { - // Coordinates with nan in S2 are unpredictable; censor to EMPTY. Empty - // points coming from WKB are always nan, nan. - if (!std::isnan(coord.x) && !std::isnan(coord.y)) { - points.push_back(S2LatLng::FromDegrees(coord.y, coord.x).Normalized().ToPoint()); - } - } - - std::unique_ptr build() { - return absl::make_unique(std::move(this->points)); - } - - private: - std::vector points; - }; - -private: - std::vector points; -}; - -#endif diff --git a/src/polygon-geography.h b/src/polygon-geography.h deleted file mode 100644 index db287718..00000000 --- a/src/polygon-geography.h +++ /dev/null @@ -1,334 +0,0 @@ - -#ifndef POLYGON_GEOGRAPHY_H -#define POLYGON_GEOGRAPHY_H - -#include "wk/reader.hpp" - -#include "geography.h" -#include "point-geography.h" -#include "polyline-geography.h" - -// This class handles polygons (POLYGON and MULTIPOLYGON) -// This is similar to an S2PolygonLayer -class PolygonGeography: public Geography { -public: - PolygonGeography() {} - PolygonGeography(std::unique_ptr polygon): - polygon(std::move(polygon)) {} - - Geography::Type GeographyType() { - return Geography::Type::GEOGRAPHY_POLYGON; - } - - bool FindValidationError(S2Error* error) { - return this->polygon->FindValidationError(error); - } - - const S2Polygon* Polygon() { - return this->polygon.get(); - } - - bool IsCollection() { - return this->outerLoopIndices().size() > 1; - } - - int Dimension() { - return 2; - } - - int NumPoints() { - return this->polygon->num_vertices(); - } - - bool IsEmpty() { - return this->polygon->is_empty(); - } - - double Area() { - return this->polygon->GetArea(); - } - - double Length() { - return 0; - } - - double Perimeter() { - std::unique_ptr boundary = this->Boundary(); - return boundary->Length(); - } - - double X() { - Rcpp::stop("Can't compute X value of a non-point geography"); - } - - double Y() { - Rcpp::stop("Can't compute Y value of a non-point geography"); - } - - S2Point Centroid() { - return this->polygon->GetCentroid(); - } - - S2Cap GetCapBound() { - return this->polygon->GetCapBound(); - } - - S2LatLngRect GetRectBound() { - return this->polygon->GetRectBound(); - } - - std::unique_ptr Boundary() { - PolylineGeography::Builder builder; - std::vector> flatIndices = this->flatLoopIndices(); - - // export multilinestring - WKGeometryMeta meta(WKGeometryType::MultiLineString, false, false, false); - meta.hasSize = true; - meta.size = this->polygon->num_loops(); - - builder.nextGeometryStart(meta, WKReader::PART_ID_NONE); - int loopId = 0; - for (size_t i = 0; i < flatIndices.size(); i++) { - this->exportLoops(&builder, meta, flatIndices[i], loopId); - loopId += flatIndices[i].size(); - } - builder.nextGeometryEnd(meta, WKReader::PART_ID_NONE); - - return builder.build(); - } - - std::vector BuildShapeIndex(MutableS2ShapeIndex* index) { - std::vector shapeIds(1); - std::unique_ptr shape = absl::make_unique(); - shape->Init(this->polygon.get()); - shapeIds[0] = index->Add(std::move(shape)); - return shapeIds; - } - - void Export(WKGeometryHandler* handler, uint32_t partId) { - std::vector> flatIndices = this->flatLoopIndices(); - - if (flatIndices.size() > 1) { - // export multipolygon - WKGeometryMeta meta(WKGeometryType::MultiPolygon, false, false, false); - meta.hasSize = true; - meta.size = flatIndices.size(); - - WKGeometryMeta childMeta(WKGeometryType::Polygon, false, false, false); - childMeta.hasSize = true; - - handler->nextGeometryStart(meta, partId); - for (size_t i = 0; i < flatIndices.size(); i++) { - childMeta.size = flatIndices[i].size(); - handler->nextGeometryStart(childMeta, i); - this->exportLoops(handler, childMeta, flatIndices[i]); - handler->nextGeometryEnd(childMeta, i); - } - - handler->nextGeometryEnd(meta, partId); - - } else if (flatIndices.size() > 0) { - // export polygon - WKGeometryMeta meta(WKGeometryType::Polygon, false, false, false); - meta.hasSize = true; - meta.size = flatIndices[0].size(); - handler->nextGeometryStart(meta, partId); - this->exportLoops(handler, meta, flatIndices[0]); - handler->nextGeometryEnd(meta, partId); - - } else { - // export empty polygon - WKGeometryMeta meta(WKGeometryType::Polygon, false, false, false); - meta.hasSize = true; - meta.size = 0; - handler->nextGeometryStart(meta, partId); - handler->nextGeometryEnd(meta, partId); - } - } - - class Builder: public GeographyBuilder { - public: - Builder(bool oriented, bool check): - oriented(oriented), check(check) {} - - void nextLinearRingStart(const WKGeometryMeta& meta, uint32_t size, uint32_t ringId) { - // skip the last vertex (WKB rings are theoretically closed) - if (size > 0) { - this->vertices = std::vector(size - 1); - } else { - this->vertices = std::vector(); - } - } - - void nextCoordinate(const WKGeometryMeta& meta, const WKCoord& coord, uint32_t coordId) { - if (coordId < this->vertices.size()) { - vertices[coordId] = S2LatLng::FromDegrees(coord.y, coord.x).Normalized().ToPoint(); - } - } - - void nextLinearRingEnd(const WKGeometryMeta& meta, uint32_t size, uint32_t ringId) { - std::unique_ptr loop = absl::make_unique(); - loop->set_s2debug_override(S2Debug::DISABLE); - loop->Init(vertices); - - if (!oriented) { - loop->Normalize(); - } - - if (this->check && !loop->IsValid()) { - std::stringstream err; - err << "Loop " << (this->loops.size()) << " is not valid: "; - S2Error error; - loop->FindValidationError(&error); - err << error.text(); - throw WKParseException(err.str()); - } - - this->loops.push_back(std::move(loop)); - } - - std::unique_ptr build() { - std::unique_ptr polygon = absl::make_unique(); - polygon->set_s2debug_override(S2Debug::DISABLE); - if (this->loops.size() > 0 && oriented) { - polygon->InitOriented(std::move(this->loops)); - } else if (this->loops.size() > 0) { - polygon->InitNested(std::move(this->loops)); - } - - // make sure polygon is valid - if (this->check && !polygon->IsValid()) { - S2Error error; - polygon->FindValidationError(&error); - throw WKParseException(error.text()); - } - - return absl::make_unique(std::move(polygon)); - } - - private: - bool oriented; - bool check; - std::vector vertices; - std::vector> loops; - }; - -private: - std::unique_ptr polygon; - - // Calculate which loops in the polygon are outer loops (loop->depth() == 0) - std::vector outerLoopIndices() { - std::vector indices; - for (int i = 0; i < this->polygon->num_loops(); i++) { - if (this->polygon->loop(i)->depth() == 0) { - indices.push_back(i); - } - } - - return indices; - } - - // Calculate the arrangement of loops in the form of a multipolygon - // (list(list(shell, !!! holes))) - std::vector> flatLoopIndices() { - std::vector outerLoops = this->outerLoopIndices(); - - std::vector> flatIndices(outerLoops.size()); - for (size_t i = 0; i < outerLoops.size(); i++) { - int k = outerLoops[i]; - flatIndices[i] = std::vector(); - - // the first loop here is the shell (depth == 0) - flatIndices[i].push_back(k); - - // loops in the S2Polygon are arranged such that child loops are - // directly after the outer loop, so add all loop indices before - // the next parent loop (or end of polygon). This is similar to - // S2Polygon::GetLastDescendant() but is slightly easier to understand. - while (++k < this->polygon->num_loops() && this->polygon->loop(k)->depth() > 0) { - flatIndices[i].push_back(k); - } - } - - return flatIndices; - } - - void exportLoops(WKGeometryHandler* handler, WKGeometryMeta meta, - const std::vector& loopIndices, int loopIdOffset = 0) { - S2LatLng point; - - for (size_t i = 0; i < loopIndices.size(); i++) { - int loopId = loopIndices[i]; - S2Loop* loop = this->polygon->loop(loopId); - if (loop->num_vertices() == 0) { - continue; - } - - // this is a slightly ugly way to make it possible to export either the - // boundaries or the loops using the same code - WKGeometryMeta childMeta(WKGeometryType::LineString, false, false, false); - childMeta.hasSize = true; - childMeta.size = loop->num_vertices() + 1; - - WKGeometryMeta coordMeta; - - if (meta.geometryType == WKGeometryType::Polygon) { - handler->nextLinearRingStart(meta, loop->num_vertices() + 1, i + loopIdOffset); - coordMeta = meta; - } else if (meta.geometryType == WKGeometryType::MultiLineString) { - handler->nextGeometryStart(childMeta, i + loopIdOffset); - coordMeta = childMeta; - } else { - std::stringstream err; - err << "Can't export S2Loop with parent geometry type " << meta.geometryType; - Rcpp::stop(err.str()); - } - - if ((loop->depth() % 2) == 0) { - // if this is the first ring, use the internal vertex order - for (int j = 0; j < loop->num_vertices(); j++) { - point = S2LatLng(loop->vertex(j)); - handler->nextCoordinate( - coordMeta, - WKCoord::xy(point.lng().degrees(), point.lat().degrees()), - j - ); - } - - // close the loop! - point = S2LatLng(loop->vertex(0)); - handler->nextCoordinate( - coordMeta, - WKCoord::xy(point.lng().degrees(), point.lat().degrees()), - loop->num_vertices() - ); - } else { - // if an interior ring, reverse the vertex order - for (int j = 0; j < loop->num_vertices(); j++) { - point = S2LatLng(loop->vertex(loop->num_vertices() - 1 - j)); - handler->nextCoordinate( - coordMeta, - WKCoord::xy(point.lng().degrees(), point.lat().degrees()), - j - ); - } - - // close the loop! - point = S2LatLng(loop->vertex(loop->num_vertices() - 1)); - handler->nextCoordinate( - coordMeta, - WKCoord::xy(point.lng().degrees(), point.lat().degrees()), - loop->num_vertices() - ); - } - - if (meta.geometryType == WKGeometryType::Polygon) { - handler->nextLinearRingEnd(meta, loop->num_vertices() + 1, i + loopIdOffset); - } else if (meta.geometryType == WKGeometryType::MultiLineString) { - handler->nextGeometryEnd(childMeta, i + loopIdOffset); - } - } - } -}; - -#endif diff --git a/src/polyline-geography.h b/src/polyline-geography.h deleted file mode 100644 index 32ffef12..00000000 --- a/src/polyline-geography.h +++ /dev/null @@ -1,215 +0,0 @@ - -#ifndef POLYLINE_GEOGRAPHY_H -#define POLYLINE_GEOGRAPHY_H - -#include "s2/s2latlng_rect.h" - -#include "geography.h" - -// This class handles (vectors of) polylines (LINESTRING and MULTILINESTRING) -// This is similar to an S2PolylineVectorLayer -class PolylineGeography: public Geography { -public: - PolylineGeography(): polylines(0) {} - PolylineGeography(std::vector> polylines): - polylines(std::move(polylines)) {} - - Geography::Type GeographyType() { - return Geography::Type::GEOGRAPHY_POLYLINE; - } - - bool FindValidationError(S2Error* error) { - bool result; - error->Clear(); - for (size_t i = 0; i < this->polylines.size(); i++) { - result = this->polylines[i]->FindValidationError(error); - if (result) { - return result; - } - } - - return false; - } - - const std::vector>* Polyline() { - return &(this->polylines); - } - - bool IsCollection() { - return this->polylines.size() > 1; - } - - int Dimension() { - return 1; - } - - int NumPoints() { - int numPoints = 0; - for (size_t i = 0; i < this->polylines.size(); i++) { - numPoints += this->polylines[i]->num_vertices(); - } - - return numPoints; - } - - bool IsEmpty() { - for (size_t i = 0; i < this->polylines.size(); i++) { - if (this->polylines[i]->num_vertices() > 0) { - return false; - } - } - - return true; - } - - double Area() { - return 0; - } - - double Length() { - double length = 0; - for (size_t i = 0; i < this->polylines.size(); i++) { - length += this->polylines[i]->GetLength().radians(); - } - - return length; - } - - double Perimeter() { - return 0; - } - - double X() { - Rcpp::stop("Can't compute X value of a non-point geography"); - } - - double Y() { - Rcpp::stop("Can't compute Y value of a non-point geography"); - } - - S2Point Centroid() { - S2Point output(0, 0, 0); - for (size_t i = 0; i < this->polylines.size(); i++) { - output += this->polylines[i]->GetCentroid(); - } - - return output; - } - - S2LatLngRect GetRectBound() { - S2LatLngRect rect; - if (this->polylines.size()) - rect = this->polylines[0]->GetRectBound(); - for (size_t i = 1; i < this->polylines.size(); i++) { - rect.Union(this->polylines[i]->GetRectBound()); // depends on order - } - return rect; - } - - std::unique_ptr Boundary() { - std::vector endpoints; - for (size_t i = 0; i < this->polylines.size(); i++) { - if (this->polylines[i]->num_vertices() >= 2) { - endpoints.push_back(this->polylines[i]->vertex(0)); - endpoints.push_back(this->polylines[i]->vertex(1)); - } - } - - return absl::make_unique(endpoints); - } - - std::vector BuildShapeIndex(MutableS2ShapeIndex* index) { - std::vector shapeIds(this->polylines.size()); - for (size_t i = 0; i < this->polylines.size(); i++) { - std::unique_ptr shape = absl::make_unique(); - shape->Init(this->polylines[i].get()); - shapeIds[i] = index->Add(std::move(shape)); - } - return shapeIds; - } - - void Export(WKGeometryHandler* handler, uint32_t partId) { - S2LatLng point; - - if (this->polylines.size() > 1) { - // export multilinestring - WKGeometryMeta meta(WKGeometryType::MultiLineString, false, false, false); - meta.hasSize = true; - meta.size = this->polylines.size(); - - handler->nextGeometryStart(meta, partId); - - for (size_t i = 0; i < this->polylines.size(); i++) { - WKGeometryMeta childMeta(WKGeometryType::LineString, false, false, false); - childMeta.hasSize = true; - childMeta.size = this->polylines[i]->num_vertices(); - - handler->nextGeometryStart(childMeta, i); - - for (size_t j = 0; j < childMeta.size; j++) { - point = S2LatLng(this->polylines[i]->vertex(j)); - handler->nextCoordinate(meta, WKCoord::xy(point.lng().degrees(), point.lat().degrees()), j); - } - - handler->nextGeometryEnd(childMeta, i); - } - - handler->nextGeometryEnd(meta, partId); - - } else if (this->polylines.size() > 0) { - // export linestring - WKGeometryMeta meta(WKGeometryType::LineString, false, false, false); - meta.hasSize = true; - meta.size = this->polylines[0]->num_vertices(); - - handler->nextGeometryStart(meta, partId); - - for (size_t i = 0; i < meta.size; i++) { - point = S2LatLng(this->polylines[0]->vertex(i)); - handler->nextCoordinate(meta, WKCoord::xy(point.lng().degrees(), point.lat().degrees()), i); - } - - handler->nextGeometryEnd(meta, partId); - - } else { - // export empty linestring - WKGeometryMeta meta(WKGeometryType::LineString, false, false, false); - meta.hasSize = true; - meta.size = this->polylines[0]->num_vertices(); - handler->nextGeometryStart(meta, partId); - handler->nextGeometryEnd(meta, partId); - } - } - - class Builder: public GeographyBuilder { - public: - void nextGeometryStart(const WKGeometryMeta& meta, uint32_t partId) { - if (meta.geometryType == WKGeometryType::LineString) { - points = std::vector(meta.size); - } - } - - void nextCoordinate(const WKGeometryMeta& meta, const WKCoord& coord, uint32_t coordId) { - points[coordId] = S2LatLng::FromDegrees(coord.y, coord.x).Normalized().ToPoint(); - } - - void nextGeometryEnd(const WKGeometryMeta& meta, uint32_t partId) { - if (meta.geometryType == WKGeometryType::LineString) { - polylines.push_back(absl::make_unique(std::move(points))); - } - } - - std::unique_ptr build() { - return absl::make_unique(std::move(this->polylines)); - } - - private: - std::vector points; - std::vector> polylines; - }; - -private: - std::vector> polylines; -}; - -#endif diff --git a/src/s2-accessors.cpp b/src/s2-accessors.cpp index 0cb9d452..c732f94d 100644 --- a/src/s2-accessors.cpp +++ b/src/s2-accessors.cpp @@ -1,15 +1,15 @@ #include "geography-operator.h" -#include "s2/s2closest_edge_query.h" -#include "s2/s2furthest_edge_query.h" #include using namespace Rcpp; +#include "s2-geography/s2-geography.hpp" + // [[Rcpp::export]] LogicalVector cpp_s2_is_collection(List geog) { class Op: public UnaryGeographyOperator { int processFeature(XPtr feature, R_xlen_t i) { - return feature->IsCollection(); + return s2geography::s2_is_collection(feature->Geog()); } }; @@ -21,7 +21,7 @@ LogicalVector cpp_s2_is_collection(List geog) { LogicalVector cpp_s2_is_valid(List geog) { class Op: public UnaryGeographyOperator { int processFeature(XPtr feature, R_xlen_t i) { - return !(feature->FindValidationError(&(this->error))); + return !s2geography::s2_find_validation_error(feature->Geog(), &error); } S2Error error; @@ -35,7 +35,7 @@ LogicalVector cpp_s2_is_valid(List geog) { CharacterVector cpp_s2_is_valid_reason(List geog) { class Op: public UnaryGeographyOperator { String processFeature(XPtr feature, R_xlen_t i) { - if (feature->FindValidationError(&(this->error))) { + if (s2geography::s2_find_validation_error(feature->Geog(), &error)) { return this->error.text(); } else { return NA_STRING; @@ -53,7 +53,7 @@ CharacterVector cpp_s2_is_valid_reason(List geog) { IntegerVector cpp_s2_dimension(List geog) { class Op: public UnaryGeographyOperator { int processFeature(XPtr feature, R_xlen_t i) { - return feature->Dimension(); + return s2geography::s2_dimension(feature->Geog()); } }; @@ -65,7 +65,7 @@ IntegerVector cpp_s2_dimension(List geog) { IntegerVector cpp_s2_num_points(List geog) { class Op: public UnaryGeographyOperator { int processFeature(XPtr feature, R_xlen_t i) { - return feature->NumPoints(); + return s2geography::s2_num_points(feature->Geog()); } }; @@ -77,7 +77,7 @@ IntegerVector cpp_s2_num_points(List geog) { LogicalVector cpp_s2_is_empty(List geog) { class Op: public UnaryGeographyOperator { int processFeature(XPtr feature, R_xlen_t i) { - return feature->IsEmpty(); + return s2geography::s2_is_empty(feature->Geog()); } }; @@ -89,7 +89,7 @@ LogicalVector cpp_s2_is_empty(List geog) { NumericVector cpp_s2_area(List geog) { class Op: public UnaryGeographyOperator { double processFeature(XPtr feature, R_xlen_t i) { - return feature->Area(); + return s2geography::s2_area(feature->Geog()); } }; @@ -101,7 +101,7 @@ NumericVector cpp_s2_area(List geog) { NumericVector cpp_s2_length(List geog) { class Op: public UnaryGeographyOperator { double processFeature(XPtr feature, R_xlen_t i) { - return feature->Length(); + return s2geography::s2_length(feature->Geog()); } }; @@ -113,7 +113,7 @@ NumericVector cpp_s2_length(List geog) { NumericVector cpp_s2_perimeter(List geog) { class Op: public UnaryGeographyOperator { double processFeature(XPtr feature, R_xlen_t i) { - return feature->Perimeter(); + return s2geography::s2_perimeter(feature->Geog()); } }; @@ -125,7 +125,11 @@ NumericVector cpp_s2_perimeter(List geog) { NumericVector cpp_s2_x(List geog) { class Op: public UnaryGeographyOperator { double processFeature(XPtr feature, R_xlen_t i) { - return feature->X(); + if (s2geography::s2_dimension(feature->Geog()) != 0) { + Rcpp::stop("Can't compute X value of a non-point geography"); + } + + return s2geography::s2_x(feature->Geog()); } }; @@ -137,7 +141,11 @@ NumericVector cpp_s2_x(List geog) { NumericVector cpp_s2_y(List geog) { class Op: public UnaryGeographyOperator { double processFeature(XPtr feature, R_xlen_t i) { - return feature->Y(); + if (s2geography::s2_dimension(feature->Geog()) != 0) { + Rcpp::stop("Can't compute Y value of a non-point geography"); + } + + return s2geography::s2_y(feature->Geog()); } }; @@ -151,27 +159,7 @@ NumericVector cpp_s2_project_normalized(List geog1, List geog2) { double processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - if (feature1->IsCollection() || feature2->IsCollection()) { - throw GeographyOperatorException("`x` and `y` must both be simple geographies"); - } - - if (feature1->IsEmpty() || feature2->IsEmpty()) { - return NA_REAL; - } - - if (feature1->GeographyType() == Geography::Type::GEOGRAPHY_POLYLINE) { - if (feature2->GeographyType() == Geography::Type::GEOGRAPHY_POINT) { - S2Point point = feature2->Point()->at(0); - int next_vertex; - S2Point point_on_line = feature1->Polyline()->at(0)->Project(point, &next_vertex); - return feature1->Polyline()->at(0)->UnInterpolate(point_on_line, next_vertex); - } else { - throw GeographyOperatorException("`y` must be a point geography"); - } - } else { - throw GeographyOperatorException("`x` must be a polyline geography"); - } - return NA_REAL; + return s2geography::s2_project_normalized(feature1->Geog(), feature2->Geog()); } }; @@ -186,13 +174,7 @@ NumericVector cpp_s2_distance(List geog1, List geog2) { double processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - S2ClosestEdgeQuery query(feature1->ShapeIndex()); - S2ClosestEdgeQuery::ShapeIndexTarget target(feature2->ShapeIndex()); - - const auto& result = query.FindClosestEdge(&target); - - S1ChordAngle angle = result.distance(); - double distance = angle.ToAngle().radians(); + double distance = s2geography::s2_distance(feature1->Index(), feature2->Index()); if (distance == R_PosInf) { return NA_REAL; @@ -213,13 +195,7 @@ NumericVector cpp_s2_max_distance(List geog1, List geog2) { double processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - S2FurthestEdgeQuery query(feature1->ShapeIndex()); - S2FurthestEdgeQuery::ShapeIndexTarget target(feature2->ShapeIndex()); - - const auto& result = query.FindFurthestEdge(&target); - - S1ChordAngle angle = result.distance(); - double distance = angle.ToAngle().radians(); + double distance = s2geography::s2_max_distance(feature1->Index(), feature2->Index()); // returns -1 if one of the indexes is empty // NA is more consistent with the BigQuery diff --git a/src/s2-bounds.cpp b/src/s2-bounds.cpp index 5e121982..e8403ad5 100644 --- a/src/s2-bounds.cpp +++ b/src/s2-bounds.cpp @@ -4,10 +4,6 @@ #include "s2-options.h" #include "geography-operator.h" -#include "point-geography.h" -#include "polyline-geography.h" -#include "polygon-geography.h" -#include "geography-collection.h" #include using namespace Rcpp; @@ -24,7 +20,7 @@ DataFrame cpp_s2_bounds_cap(List geog) { lat[i] = lng[i] = angle[i] = NA_REAL; } else { Rcpp::XPtr feature(item); - S2Cap cap = feature->GetCapBound(); + S2Cap cap = feature->Geog().Region()->GetCapBound(); S2LatLng center(cap.center()); lng[i] = center.lng().degrees(); lat[i] = center.lat().degrees(); @@ -51,7 +47,7 @@ DataFrame cpp_s2_bounds_rect(List geog) { lng_lo[i] = lat_lo[i] = lng_hi[i] = lat_hi[i] = NA_REAL; } else { Rcpp::XPtr feature(item); - S2LatLngRect rect = feature->GetRectBound(); + S2LatLngRect rect = feature->Geog().Region()->GetRectBound(); lng_lo[i] = rect.lng_lo().degrees(); lat_lo[i] = rect.lat_lo().degrees(); lng_hi[i] = rect.lng_hi().degrees(); diff --git a/src/s2-cell-union.cpp b/src/s2-cell-union.cpp index 38bc6840..a3b4adf5 100644 --- a/src/s2-cell-union.cpp +++ b/src/s2-cell-union.cpp @@ -8,9 +8,6 @@ #include "s2/s2region_union.h" #include "geography-operator.h" -#include "point-geography.h" -#include "polyline-geography.h" -#include "polygon-geography.h" #include using namespace Rcpp; @@ -308,7 +305,7 @@ List cpp_s2_geography_from_cell_union(List cellUnionVector) { SEXP processCell(S2CellUnion& cellUnion, R_xlen_t i) { std::unique_ptr polygon = absl::make_unique(); polygon->InitToCellUnionBorder(cellUnion); - return XPtr(new PolygonGeography(std::move(polygon))); + return Geography::MakeXPtr(Geography::MakePolygon(std::move(polygon))); } }; @@ -331,7 +328,7 @@ List cpp_s2_covering_cell_ids(List geog, int min_level, int max_level, SEXP processFeature(XPtr feature, R_xlen_t i) { S2ShapeIndexBufferedRegion region; - region.Init(feature->ShapeIndex(), S1ChordAngle::Radians(this->distance[i])); + region.Init(&feature->Index().ShapeIndex(), S1ChordAngle::Radians(this->distance[i])); S2CellUnion cellUnion; if (interior) { @@ -366,7 +363,6 @@ List cpp_s2_covering_cell_ids_agg(List geog, int min_level, int max_level, S1ChordAngle bufferAngle = S1ChordAngle::Radians(buffer); S2RegionUnion regionUnion; - SEXP item; for (R_xlen_t i = 0; i < geog.size(); i++) { item = geog[i]; @@ -379,7 +375,7 @@ List cpp_s2_covering_cell_ids_agg(List geog, int min_level, int max_level, if (item != R_NilValue) { Rcpp::XPtr feature(item); auto region = absl::make_unique(); - region->Init(feature->ShapeIndex(), bufferAngle); + region->Init(&feature->Index().ShapeIndex(), bufferAngle); regionUnion.Add(std::move(region)); } } diff --git a/src/s2-cell.cpp b/src/s2-cell.cpp index 12251166..f2257485 100644 --- a/src/s2-cell.cpp +++ b/src/s2-cell.cpp @@ -9,9 +9,7 @@ #include "s2/s2cell.h" #include "s2/s2latlng.h" -#include "point-geography.h" -#include "polyline-geography.h" -#include "polygon-geography.h" +#include "geography.h" #include using namespace Rcpp; @@ -356,7 +354,7 @@ List cpp_s2_cell_center(NumericVector cellIdVector) { class Op: public UnaryS2CellOperator { SEXP processCell(S2CellId cellId, R_xlen_t i) { if (cellId.is_valid()) { - return XPtr(new PointGeography(cellId.ToPoint())); + return Geography::MakeXPtr(Geography::MakePoint(cellId.ToPoint())); } else { return R_NilValue; } @@ -374,7 +372,8 @@ List cpp_s2_cell_polygon(NumericVector cellIdVector) { class Op: public UnaryS2CellOperator { SEXP processCell(S2CellId cellId, R_xlen_t i) { if (cellId.is_valid()) { - return XPtr(new PolygonGeography(absl::make_unique(S2Cell(cellId)))); + auto poly = absl::make_unique(S2Cell(cellId)); + return Geography::MakeXPtr(Geography::MakePolygon(std::move(poly))); } else { return R_NilValue; } @@ -392,7 +391,7 @@ List cpp_s2_cell_vertex(NumericVector cellIdVector, IntegerVector k) { class Op: public UnaryS2CellOperator { SEXP processCell(S2CellId cellId, R_xlen_t i) { if (cellId.is_valid() && (this->k[i] >= 0)) { - return XPtr(new PointGeography(S2Cell(cellId).GetVertex(this->k[i]))); + return Geography::MakeXPtr(Geography::MakePoint(S2Cell(cellId).GetVertex(this->k[i]))); } else { return R_NilValue; } diff --git a/src/s2-constructors-formatters.cpp b/src/s2-constructors-formatters.cpp index 7741e4d6..9aea1e1d 100644 --- a/src/s2-constructors-formatters.cpp +++ b/src/s2-constructors-formatters.cpp @@ -1,82 +1,630 @@ -#include -#include "wk/rcpp-translate.hpp" -#include "wk/rcpp-coord-reader.hpp" +#define R_NO_REMAP +#include +#include -#include "wk-geography.h" +#include "wk-v1.h" +#include "s2-geography/s2-geography.hpp" +#include "geography.h" -using namespace Rcpp; -// [[Rcpp::export]] -List cpp_s2_geog_point(NumericVector x, NumericVector y) { - NumericVector z(x.size()); - z.fill(NA_REAL); - NumericVector m(x.size()); - m.fill(NA_REAL); +#define CPP_START \ + char cpp_exception_error[8096]; \ + memset(cpp_exception_error, 0, 8096); \ + try { - WKRcppPointCoordProvider provider(x, y, z, m); - WKRcppPointCoordReader reader(provider); +#define CPP_END \ + } catch (std::exception& e) { \ + strncpy(cpp_exception_error, e.what(), 8096 - 1); \ + } \ + Rf_error("%s", cpp_exception_error); \ + return R_NilValue; - WKGeographyWriter writer(provider.nFeatures()); - reader.setHandler(&writer); - while (reader.hasNextFeature()) { - checkUserInterrupt(); - reader.iterateFeature(); +// The other versions of CPP_START and CPP_END stack-allocate the +// error message buffer, which takes a non-trivial amount of time +// when done at this scale (at worst 4 times per coordinate). By +// keeping the buffer in the handler_data struct, we can call C++ +// from every handler method without measurable overhead. +#define WK_METHOD_CPP_START \ + try { + +#define WK_METHOD_CPP_END \ + } catch (std::exception& e) { \ + strncpy(data->cpp_exception_error, e.what(), 8096 - 1); \ + } \ + Rf_error("%s", data->cpp_exception_error); \ + return R_NilValue; + +#define WK_METHOD_CPP_END_INT \ + } catch (std::exception& e) { \ + strncpy(data->cpp_exception_error, e.what(), 8096 - 1); \ + } \ + Rf_error("%s", data->cpp_exception_error); \ + return WK_ABORT; + + +typedef struct { + s2geography::VectorConstructor* builder; + SEXP result; + R_xlen_t feat_id; + int coord_size; + char cpp_exception_error[8096]; +} builder_handler_t; + + +// TODO: Both of these allocate in a way that could longjmp and possibly leak memory +static inline void builder_result_append(builder_handler_t* data, SEXP value) { + R_xlen_t current_size = Rf_xlength(data->result); + if (data->feat_id >= current_size) { + SEXP new_result = PROTECT(Rf_allocVector(VECSXP, current_size * 2 + 1)); + for (R_xlen_t i = 0; i < current_size; i++) { + SET_VECTOR_ELT(new_result, i, VECTOR_ELT(data->result, i)); + } + R_ReleaseObject(data->result); + data->result = new_result; + R_PreserveObject(data->result); + UNPROTECT(1); + } + + SET_VECTOR_ELT(data->result, data->feat_id, value); + data->feat_id++; +} + +static inline void builder_result_finalize(builder_handler_t* data) { + R_xlen_t current_size = Rf_xlength(data->result); + if (data->feat_id != current_size) { + SEXP new_result = PROTECT(Rf_allocVector(VECSXP, data->feat_id)); + for (R_xlen_t i = 0; i < data->feat_id; i++) { + SET_VECTOR_ELT(new_result, i, VECTOR_ELT(data->result, i)); + } + R_ReleaseObject(data->result); + data->result = new_result; + R_PreserveObject(data->result); + UNPROTECT(1); + } +} + +int builder_vector_start(const wk_vector_meta_t* meta, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + + if (data->result != R_NilValue) { + Rf_error("Destination vector was already allocated"); // # nocov + } + + if (meta->size == WK_VECTOR_SIZE_UNKNOWN) { + data->result = PROTECT(Rf_allocVector(VECSXP, 1024)); + } else { + data->result = PROTECT(Rf_allocVector(VECSXP, meta->size)); } + R_PreserveObject(data->result); + UNPROTECT(1); + + data->feat_id = 0; + + return WK_CONTINUE; +} + +SEXP builder_vector_end(const wk_vector_meta_t* meta, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + builder_result_finalize(data); + SEXP cls = PROTECT(Rf_allocVector(STRSXP, 2)); + SET_STRING_ELT(cls, 0, Rf_mkChar("s2_geography")); + SET_STRING_ELT(cls, 1, Rf_mkChar("s2_xptr")); + Rf_setAttrib(data->result, R_ClassSymbol, cls); + UNPROTECT(1); + return data->result; +} - return writer.output; +int builder_feature_start(const wk_vector_meta_t* meta, R_xlen_t feat_id, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + WK_METHOD_CPP_START + data->builder->start_feature(); + return WK_CONTINUE; + WK_METHOD_CPP_END_INT } -// [[Rcpp::export]] -List cpp_s2_make_line(NumericVector x, NumericVector y, IntegerVector featureId) { - NumericVector z(x.size()); - z.fill(NA_REAL); - NumericVector m(x.size()); - m.fill(NA_REAL); +int builder_feature_null(void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + builder_result_append(data, R_NilValue); + return WK_ABORT_FEATURE; +} - WKRcppLinestringCoordProvider provider(x, y, z, m, featureId); - WKRcppLinestringCoordReader reader(provider); +int builder_feature_end(const wk_vector_meta_t* meta, R_xlen_t feat_id, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + WK_METHOD_CPP_START + std::unique_ptr feat = data->builder->finish_feature(); + builder_result_append(data, Geography::MakeXPtr(std::move(feat))); + return WK_CONTINUE; + WK_METHOD_CPP_END_INT +} - WKGeographyWriter writer(provider.nFeatures()); - reader.setHandler(&writer); +int builder_geometry_start(const wk_meta_t* meta, uint32_t part_id, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + WK_METHOD_CPP_START - while (reader.hasNextFeature()) { - checkUserInterrupt(); - reader.iterateFeature(); + auto geometry_type = static_cast(meta->geometry_type); + + int32_t size; + if (meta->size == WK_SIZE_UNKNOWN) { + size = -1; + } else { + size = meta->size; } - return writer.output; + if (meta->flags & WK_FLAG_HAS_Z && meta->flags & WK_FLAG_HAS_M) { + data->coord_size = 4; + } else if (meta->flags & WK_FLAG_HAS_Z) { + data->coord_size = 3; + } else if (meta->flags & WK_FLAG_HAS_M) { + data->coord_size = 3; + } else { + data->coord_size = 2; + } + + data->builder->geom_start(geometry_type, size); + return WK_CONTINUE; + WK_METHOD_CPP_END_INT +} + +int builder_geometry_end(const wk_meta_t* meta, uint32_t part_id, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + WK_METHOD_CPP_START + data->builder->geom_end(); + return WK_CONTINUE; + WK_METHOD_CPP_END_INT } -// [[Rcpp::export]] -List cpp_s2_make_polygon(NumericVector x, NumericVector y, - IntegerVector featureId, IntegerVector ringId, - bool oriented, bool check) { - NumericVector z(x.size()); - z.fill(NA_REAL); - NumericVector m(x.size()); - m.fill(NA_REAL); +int builder_ring_start(const wk_meta_t* meta, uint32_t size, uint32_t ring_id, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + WK_METHOD_CPP_START + + if (size == WK_SIZE_UNKNOWN) { + data->builder->ring_start(-1); + } else { + data->builder->ring_start(size); + } - WKRcppPolygonCoordProvider provider(x, y, z, m, featureId, ringId); - WKRcppPolygonCoordReader reader(provider); + return WK_CONTINUE; + WK_METHOD_CPP_END_INT +} - WKGeographyWriter writer(provider.nFeatures()); - writer.setOriented(oriented); - writer.setCheck(check); +int builder_ring_end(const wk_meta_t* meta, uint32_t size, uint32_t ring_id, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + WK_METHOD_CPP_START + data->builder->ring_end(); + return WK_CONTINUE; + WK_METHOD_CPP_END_INT +} - reader.setHandler(&writer); +int builder_coord(const wk_meta_t* meta, const double* coord, uint32_t coord_id, void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + WK_METHOD_CPP_START + data->builder->coords(coord, 1, data->coord_size); + return WK_CONTINUE; + WK_METHOD_CPP_END_INT +} - while (reader.hasNextFeature()) { - checkUserInterrupt(); - reader.iterateFeature(); +int builder_error(const char* message, void* handler_data) { + Rf_error("%s", message); + return WK_ABORT; +} + +void builder_deinitialize(void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + if (data->result != R_NilValue) { + R_ReleaseObject(data->result); + data->result = R_NilValue; } +} - if (writer.problemId.size() > 0) { - Environment s2NS = Environment::namespace_env("s2"); - Function stopProblems = s2NS["stop_problems_create"]; - stopProblems(writer.problemId, writer.problems); +void builder_finalize(void* handler_data) { + builder_handler_t* data = (builder_handler_t*) handler_data; + if (data != nullptr) { + free(data); } +} + +void delete_vector_constructor(SEXP xptr) { + auto ptr = reinterpret_cast(R_ExternalPtrAddr(xptr)); + if (ptr != nullptr) { + delete ptr; + } +} + +extern "C" SEXP c_s2_geography_writer_new(SEXP oriented_sexp, SEXP check_sexp) { + CPP_START + + int oriented = LOGICAL(oriented_sexp)[0]; + int check = LOGICAL(check_sexp)[0]; + + s2geography::Constructor::Options options; + options.set_oriented(oriented); + options.set_check(check); + + auto builder = new s2geography::VectorConstructor(options); + SEXP builder_xptr = PROTECT(R_MakeExternalPtr(builder, R_NilValue, R_NilValue)); + R_RegisterCFinalizer(builder_xptr, &delete_vector_constructor); + + wk_handler_t* handler = wk_handler_create(); + + handler->vector_start = &builder_vector_start; + handler->vector_end = &builder_vector_end; + + handler->feature_start = &builder_feature_start; + handler->null_feature = &builder_feature_null; + handler->feature_end = &builder_feature_end; + + handler->geometry_start = &builder_geometry_start; + handler->geometry_end = &builder_geometry_end; + + handler->ring_start = &builder_ring_start; + handler->ring_end = &builder_ring_end; + + handler->coord = &builder_coord; + + handler->error = &builder_error; + + handler->deinitialize = &builder_deinitialize; + handler->finalizer = &builder_finalize; + + builder_handler_t* data = (builder_handler_t*) malloc(sizeof(builder_handler_t)); + if (data == NULL) { + wk_handler_destroy(handler); // # nocov + Rf_error("Failed to alloc handler data"); // # nocov + } + + data->coord_size = 2; + data->builder = builder; + data->result = R_NilValue; + memset(data->cpp_exception_error, 0, 8096); + + handler->handler_data = data; + + // include the builder pointer as a tag for this external pointer + // which guarnatees that it will not be garbage collected until + // this object is garbage collected + SEXP handler_xptr = wk_handler_create_xptr(handler, builder_xptr, R_NilValue); + UNPROTECT(1); + return handler_xptr; + + CPP_END +} + +#define HANDLE_OR_RETURN(expr) \ + result = expr; \ + if (result != WK_CONTINUE) return result + +#define HANDLE_CONTINUE_OR_BREAK(expr) \ + result = expr; \ + if (result == WK_ABORT_FEATURE) continue; else if (result == WK_ABORT) break + + +int handle_points(const s2geography::PointGeography& geog, wk_handler_t* handler, + uint32_t part_id = WK_PART_ID_NONE) { + int result; + + wk_meta_t meta; + WK_META_RESET(meta, WK_MULTIPOINT); + meta.size = geog.Points().size(); + + wk_meta_t meta_child; + WK_META_RESET(meta_child, WK_POINT); + meta_child.size = 1; + double coord[2]; + + if (meta.size == 0) { + meta_child.size = 0; + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, part_id, handler->handler_data)); + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, part_id, handler->handler_data)); + } else if (meta.size == 1) { + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, part_id, handler->handler_data)); + S2LatLng pt(geog.Points()[0]); + coord[0] = pt.lng().degrees(); + coord[1] = pt.lat().degrees(); + HANDLE_OR_RETURN(handler->coord(&meta_child, coord, 0, handler->handler_data)); + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, part_id, handler->handler_data)); + } else { + HANDLE_OR_RETURN(handler->geometry_start(&meta, part_id, handler->handler_data)); + + for (size_t i = 0; i < geog.Points().size(); i++) { + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, i, handler->handler_data)); + S2LatLng pt(geog.Points()[i]); + coord[0] = pt.lng().degrees(); + coord[1] = pt.lat().degrees(); + HANDLE_OR_RETURN(handler->coord(&meta_child, coord, 0, handler->handler_data)); + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, i, handler->handler_data)); + } + + HANDLE_OR_RETURN(handler->geometry_end(&meta, part_id, handler->handler_data)); + } + + return WK_CONTINUE; +} + +int handle_polylines(const s2geography::PolylineGeography& geog, wk_handler_t* handler, + uint32_t part_id = WK_PART_ID_NONE) { + int result; + + wk_meta_t meta; + WK_META_RESET(meta, WK_MULTILINESTRING); + meta.size = geog.Polylines().size(); + + wk_meta_t meta_child; + WK_META_RESET(meta_child, WK_LINESTRING); + double coord[2]; + + if (meta.size == 0) { + meta_child.size = 0; + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, part_id, handler->handler_data)); + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, part_id, handler->handler_data)); + } else if (meta.size == 1) { + const S2Polyline& poly = *geog.Polylines()[0]; + meta_child.size = poly.num_vertices(); + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, part_id, handler->handler_data)); + + for (int j = 0; j < poly.num_vertices(); j++) { + S2LatLng pt(poly.vertex(j)); + coord[0] = pt.lng().degrees(); + coord[1] = pt.lat().degrees(); + HANDLE_OR_RETURN(handler->coord(&meta_child, coord, j, handler->handler_data)); + } + + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, part_id, handler->handler_data)); + } else { + HANDLE_OR_RETURN(handler->geometry_start(&meta, part_id, handler->handler_data)); + + for (size_t i = 0; i < geog.Polylines().size(); i++) { + const S2Polyline& poly = *geog.Polylines()[i]; + meta_child.size = poly.num_vertices(); + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, i, handler->handler_data)); + + for (int j = 0; j < poly.num_vertices(); j++) { + S2LatLng pt(poly.vertex(j)); + coord[0] = pt.lng().degrees(); + coord[1] = pt.lat().degrees(); + HANDLE_OR_RETURN(handler->coord(&meta_child, coord, j, handler->handler_data)); + } + + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, i, handler->handler_data)); + } + + HANDLE_OR_RETURN(handler->geometry_end(&meta, part_id, handler->handler_data)); + } + + return WK_CONTINUE; +} + +int handle_loop_shell(const S2Loop* loop, const wk_meta_t* meta, uint32_t loop_id, wk_handler_t* handler) { + int result; + double coord[2]; + + if (loop->num_vertices() == 0) { + return handler->error("Unexpected S2Loop with 0 vertices", handler->handler_data); + } + + HANDLE_OR_RETURN(handler->ring_start(meta, loop->num_vertices() + 1, loop_id, handler->handler_data)); + + for (int i = 0; i <= loop->num_vertices(); i++) { + S2LatLng pt(loop->vertex(i)); + coord[0] = pt.lng().degrees(); + coord[1] = pt.lat().degrees(); + HANDLE_OR_RETURN(handler->coord(meta, coord, i, handler->handler_data)); + } + + HANDLE_OR_RETURN(handler->ring_end(meta, loop->num_vertices() + 1, loop_id, handler->handler_data)); + return WK_CONTINUE; +} + +int handle_loop_hole(const S2Loop* loop, const wk_meta_t* meta, uint32_t loop_id, wk_handler_t* handler) { + int result; + double coord[2]; + + if (loop->num_vertices() == 0) { + return handler->error("Unexpected S2Loop with 0 vertices", handler->handler_data); + } + + HANDLE_OR_RETURN(handler->ring_start(meta, loop->num_vertices() + 1, loop_id, handler->handler_data)); + + uint32_t coord_id = 0; + for (int i = loop->num_vertices() - 1; i >= 0; i--) { + S2LatLng pt(loop->vertex(i)); + coord[0] = pt.lng().degrees(); + coord[1] = pt.lat().degrees(); + HANDLE_OR_RETURN(handler->coord(meta, coord, coord_id, handler->handler_data)); + coord_id++; + } + + S2LatLng pt(loop->vertex(loop->num_vertices() - 1)); + coord[0] = pt.lng().degrees(); + coord[1] = pt.lat().degrees(); + HANDLE_OR_RETURN(handler->coord(meta, coord, coord_id, handler->handler_data)); + + HANDLE_OR_RETURN(handler->ring_end(meta, loop->num_vertices() + 1, loop_id, handler->handler_data)); + return WK_CONTINUE; +} + +int handle_shell(const S2Polygon& poly, const wk_meta_t* meta, int loop_start, wk_handler_t* handler) { + int result; + const S2Loop* loop0 = poly.loop(loop_start); + HANDLE_OR_RETURN(handle_loop_shell(loop0, meta, 0, handler)); + + uint32_t loop_id = 1; + for (int j = loop_start + 1; j <= poly.GetLastDescendant(loop_start); j++) { + const S2Loop* loop = poly.loop(j); + if (loop->depth() == (loop0->depth() + 1)) { + HANDLE_OR_RETURN(handle_loop_hole(loop, meta, loop_id, handler)); + loop_id++; + } + } + + return WK_CONTINUE; +} + +int handle_polygon(const s2geography::PolygonGeography& geog, wk_handler_t* handler, + uint32_t part_id = WK_PART_ID_NONE) { + const S2Polygon& poly = *geog.Polygon(); + + // find the outer shells (loop depth = 0, 2, 4, etc.) + std::vector outer_shell_loop_ids; + std::vector outer_shell_loop_sizes; + + outer_shell_loop_ids.reserve(poly.num_loops()); + for (int i = 0; i < poly.num_loops(); i++) { + if ((poly.loop(i)->depth() % 2) == 0) { + outer_shell_loop_ids.push_back(i); + } + } + + // count the number of rings in each + outer_shell_loop_sizes.reserve(outer_shell_loop_ids.size()); + for (const auto loop_start : outer_shell_loop_ids) { + const S2Loop* loop0 = poly.loop(loop_start); + int num_loops = 1; + + for (int j = loop_start + 1; j <= poly.GetLastDescendant(loop_start); j++) { + const S2Loop* loop = poly.loop(j); + num_loops += loop->depth() == (loop0->depth() + 1); + } + + outer_shell_loop_sizes.push_back(num_loops); + } + + int result; + + wk_meta_t meta; + WK_META_RESET(meta, WK_MULTIPOLYGON); + meta.size = outer_shell_loop_ids.size(); + + wk_meta_t meta_child; + WK_META_RESET(meta_child, WK_POLYGON); + + if (meta.size == 0) { + meta_child.size = 0; + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, part_id, handler->handler_data)); + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, part_id, handler->handler_data)); + } else if (meta.size == 1) { + meta_child.size = outer_shell_loop_sizes[0]; + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, part_id, handler->handler_data)); + HANDLE_OR_RETURN(handle_shell(poly, &meta_child, outer_shell_loop_ids[0], handler)); + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, part_id, handler->handler_data)); + } else { + HANDLE_OR_RETURN(handler->geometry_start(&meta, part_id, handler->handler_data)); + + for (size_t i = 0; i < outer_shell_loop_sizes.size(); i++) { + meta_child.size = outer_shell_loop_sizes[i]; + HANDLE_OR_RETURN(handler->geometry_start(&meta_child, i, handler->handler_data)); + HANDLE_OR_RETURN(handle_shell(poly, &meta_child, outer_shell_loop_ids[i], handler)); + HANDLE_OR_RETURN(handler->geometry_end(&meta_child, i, handler->handler_data)); + } + + HANDLE_OR_RETURN(handler->geometry_end(&meta, part_id, handler->handler_data)); + } + + return WK_CONTINUE; +} + +int handle_collection(const s2geography::S2GeographyCollection& geog, wk_handler_t* handler, + uint32_t part_id = WK_PART_ID_NONE) { + int result; + + wk_meta_t meta; + WK_META_RESET(meta, WK_GEOMETRYCOLLECTION); + meta.size = geog.Features().size(); + + HANDLE_OR_RETURN(handler->geometry_start(&meta, part_id, handler->handler_data)); + for (size_t i = 0; i < geog.Features().size(); i++) { + const s2geography::S2Geography* child_ptr = geog.Features()[i].get(); + + auto child_point = dynamic_cast(child_ptr); + if (child_point != nullptr) { + HANDLE_OR_RETURN(handle_points(*child_point, handler, i)); + continue; + } + + auto child_polyline = dynamic_cast(child_ptr); + if (child_polyline != nullptr) { + HANDLE_OR_RETURN(handle_polylines(*child_polyline, handler, i)); + continue; + } + + auto child_polygon = dynamic_cast(child_ptr); + if (child_polygon != nullptr) { + HANDLE_OR_RETURN(handle_polygon(*child_polygon, handler, i)); + continue; + } + + auto child_collection = dynamic_cast(child_ptr); + if (child_collection != nullptr) { + HANDLE_OR_RETURN(handle_collection(*child_collection, handler, i)); + continue; + } + + return handler->error("Unsupported S2Geography subclass", handler->handler_data); + } + HANDLE_OR_RETURN(handler->geometry_end(&meta, part_id, handler->handler_data)); + + return WK_CONTINUE; +} + +SEXP handle_geography(SEXP data, wk_handler_t* handler) { + R_xlen_t n_features = Rf_xlength(data); + + wk_vector_meta_t vector_meta; + WK_VECTOR_META_RESET(vector_meta, WK_GEOMETRY); + vector_meta.size = n_features; + vector_meta.flags |= WK_FLAG_DIMS_UNKNOWN; + + if (handler->vector_start(&vector_meta, handler->handler_data) == WK_CONTINUE) { + int result; + SEXP item; + + for (R_xlen_t i = 0; i < n_features; i++) { + item = VECTOR_ELT(data, i); + + HANDLE_CONTINUE_OR_BREAK(handler->feature_start(&vector_meta, i, handler->handler_data)); + + if (item == R_NilValue) { + HANDLE_CONTINUE_OR_BREAK(handler->null_feature(handler->handler_data)); + } else { + auto item_ptr = reinterpret_cast(R_ExternalPtrAddr(item)); + const s2geography::S2Geography* geog_ptr = &item_ptr->Geog(); + + auto child_point = dynamic_cast(geog_ptr); + if (child_point != nullptr) { + HANDLE_CONTINUE_OR_BREAK(handle_points(*child_point, handler)); + } else { + auto child_polyline = dynamic_cast(geog_ptr); + if (child_polyline != nullptr) { + HANDLE_CONTINUE_OR_BREAK(handle_polylines(*child_polyline, handler)); + } else { + auto child_polygon = dynamic_cast(geog_ptr); + if (child_polygon != nullptr) { + HANDLE_CONTINUE_OR_BREAK(handle_polygon(*child_polygon, handler)); + } else { + auto child_collection = dynamic_cast(geog_ptr); + if (child_collection != nullptr) { + HANDLE_CONTINUE_OR_BREAK(handle_collection(*child_collection, handler)); + } else { + HANDLE_CONTINUE_OR_BREAK( + handler->error("Unsupported S2Geography subclass", handler->handler_data)); + } + } + } + } + } + + if (handler->feature_end(&vector_meta, i, handler->handler_data) == WK_ABORT) { + break; + } + } + } + + SEXP result = PROTECT(handler->vector_end(&vector_meta, handler->handler_data)); + UNPROTECT(1); + return result; +} - return writer.output; +extern "C" SEXP c_s2_handle_geography(SEXP data, SEXP handler_xptr) { + return wk_handler_run_xptr(&handle_geography, data, handler_xptr); } diff --git a/src/s2-geography.cpp b/src/s2-geography.cpp index 88eb386d..2f062152 100644 --- a/src/s2-geography.cpp +++ b/src/s2-geography.cpp @@ -3,133 +3,14 @@ #include "s2/s2polyline.h" #include "s2/s2polygon.h" -#include "wk/wkb-reader.hpp" -#include "wk/wkt-reader.hpp" -#include "wk/wkb-writer.hpp" -#include "wk/wkt-writer.hpp" -#include "wk/geometry-formatter.hpp" - #include "geography.h" -#include "wk-geography.h" -#include "point-geography.h" -#include "polyline-geography.h" -#include "polygon-geography.h" -#include "geography-collection.h" #include using namespace Rcpp; - -// [[Rcpp::export]] -List s2_geography_from_wkb(List wkb, bool oriented, bool check) { - WKRawVectorListProvider provider(wkb); - WKGeographyWriter writer(wkb.size()); - writer.setOriented(oriented); - writer.setCheck(check); - - WKBReader reader(provider); - reader.setHandler(&writer); - - while (reader.hasNextFeature()) { - checkUserInterrupt(); - reader.iterateFeature(); - } - - if (writer.problemId.size() > 0) { - Environment s2NS = Environment::namespace_env("s2"); - Function stopProblems = s2NS["stop_problems_create"]; - stopProblems(writer.problemId, writer.problems); - } - - return writer.output; -} - -// [[Rcpp::export]] -List s2_geography_from_wkt(CharacterVector wkt, bool oriented, bool check) { - WKCharacterVectorProvider provider(wkt); - WKGeographyWriter writer(wkt.size()); - writer.setOriented(oriented); - writer.setCheck(check); - - WKTReader reader(provider); - reader.setHandler(&writer); - - while (reader.hasNextFeature()) { - checkUserInterrupt(); - reader.iterateFeature(); - } - - if (writer.problemId.size() > 0) { - Environment s2NS = Environment::namespace_env("s2"); - Function stopProblems = s2NS["stop_problems_create"]; - stopProblems(writer.problemId, writer.problems); - } - - return writer.output; -} - // [[Rcpp::export]] List s2_geography_full(LogicalVector x) { // create single geography with full polygon std::unique_ptr l = absl::make_unique(S2Loop::kFull()); std::unique_ptr p = absl::make_unique(std::move(l)); - Geography *pg = new PolygonGeography(std::move(p)); - List ret(1); - ret(0) = Rcpp::XPtr(pg); - return ret; -} - -// [[Rcpp::export]] -CharacterVector s2_geography_to_wkt(List s2_geography, int precision, bool trim) { - WKRcppSEXPProvider provider(s2_geography); - WKGeographyReader reader(provider); - - WKCharacterVectorExporter exporter(reader.nFeatures()); - exporter.setRoundingPrecision(precision); - exporter.setTrim(trim); - WKTWriter writer(exporter); - - reader.setHandler(&writer); - while (reader.hasNextFeature()) { - checkUserInterrupt(); - reader.iterateFeature(); - } - - return exporter.output; -} - -// [[Rcpp::export]] -List s2_geography_to_wkb(List s2_geography, int endian) { - WKRcppSEXPProvider provider(s2_geography); - WKGeographyReader reader(provider); - - WKRawVectorListExporter exporter(reader.nFeatures()); - WKBWriter writer(exporter); - writer.setEndian(endian); - - reader.setHandler(&writer); - while (reader.hasNextFeature()) { - checkUserInterrupt(); - reader.iterateFeature(); - } - - return exporter.output; -} - -// [[Rcpp::export]] -CharacterVector s2_geography_format(List s2_geography, int maxCoords, int precision, bool trim) { - WKRcppSEXPProvider provider(s2_geography); - WKGeographyReader reader(provider); - - WKCharacterVectorExporter exporter(s2_geography.size()); - exporter.setRoundingPrecision(precision); - exporter.setTrim(trim); - WKGeometryFormatter formatter(exporter, maxCoords); - - reader.setHandler(&formatter); - while (reader.hasNextFeature()) { - checkUserInterrupt(); - reader.iterateFeature(); - } - - return exporter.output; + return List::create(Geography::MakeXPtr(Geography::MakePolygon(std::move(p)))); } diff --git a/src/s2-geography/accessors-geog.cpp b/src/s2-geography/accessors-geog.cpp new file mode 100644 index 00000000..200a1546 --- /dev/null +++ b/src/s2-geography/accessors-geog.cpp @@ -0,0 +1,201 @@ + +#include "s2/s2centroids.h" + +#include "accessors-geog.hpp" +#include "geography.hpp" +#include "build.hpp" +#include "accessors.hpp" + +namespace s2geography { + +S2Point s2_centroid(const S2Geography& geog) { + S2Point centroid(0, 0, 0); + + if (geog.dimension() == 0) { + for (int i = 0; i < geog.num_shapes(); i++) { + auto shape = geog.Shape(i); + for (int j = 0; j < shape->num_edges(); j++) { + centroid += shape->edge(j).v0; + } + } + + return centroid.Normalize(); + } + + if (geog.dimension() == 1) { + for (int i = 0; i < geog.num_shapes(); i++) { + auto shape = geog.Shape(i); + for (int j = 0; j < shape->num_edges(); j++) { + S2Shape::Edge e = shape->edge(j); + centroid += S2::TrueCentroid(e.v0, e.v1); + } + } + + return centroid.Normalize(); + } + + if (geog.dimension() == 2) { + auto polygon_ptr = dynamic_cast(&geog); + if (polygon_ptr != nullptr) { + centroid = polygon_ptr->Polygon()->GetCentroid(); + } else { + std::unique_ptr built = s2_build_polygon(geog); + centroid = built->Polygon()->GetCentroid(); + } + + return centroid.Normalize(); + } + + auto collection_ptr = dynamic_cast(&geog); + if (collection_ptr == nullptr) { + throw Exception("Can't compute s2_centroid() on custom collection geography"); + } + + for (auto& feat: collection_ptr->Features()) { + centroid += s2_centroid(*feat); + } + + return centroid.Normalize(); +} + +std::unique_ptr s2_boundary(const S2Geography& geog) { + int dimension = s2_dimension(geog); + + if (dimension == 1) { + std::vector endpoints; + for (int i = 0; i < geog.num_shapes(); i++) { + auto shape = geog.Shape(i); + if (shape->dimension() < 1) { + continue; + } + + endpoints.reserve(endpoints.size() + shape->num_chains() * 2); + for (int j = 0; j < shape->num_chains(); j++) { + S2Shape::Chain chain = shape->chain(j); + if (chain.length > 0) { + endpoints.push_back(shape->edge(chain.start).v0); + endpoints.push_back(shape->edge(chain.start + chain.length - 1).v1); + } + } + } + + return absl::make_unique(std::move(endpoints)); + } + + if (dimension == 2) { + std::vector> polylines; + polylines.reserve(geog.num_shapes()); + + for (int i = 0; i < geog.num_shapes(); i++) { + auto shape = geog.Shape(i); + if (shape->dimension() != 2) { + throw Exception("Can't extract boundary from heterogeneous collection"); + } + + for (int j = 0; j < shape->num_chains(); j++) { + S2Shape::Chain chain = shape->chain(j); + if (chain.length > 0) { + std::vector points(chain.length + 1); + + points[0] = shape->edge(chain.start).v0; + for (int k = 0; k < chain.length; k++) { + points[k + 1] = shape->edge(chain.start + k).v1; + } + + auto polyline = absl::make_unique(std::move(points)); + polylines.push_back(std::move(polyline)); + } + } + } + + return absl::make_unique(std::move(polylines)); + } + + return absl::make_unique(); +} + +std::unique_ptr s2_convex_hull(const S2Geography& geog) { + S2ConvexHullAggregator agg; + agg.Add(geog); + return agg.Finalize(); +} + + +void CentroidAggregator::Add(const S2Geography& geog) { + S2Point centroid = s2_centroid(geog); + if (centroid.Norm2() > 0) { + centroid_ += centroid.Normalize(); + } +} + +void CentroidAggregator::Merge(const CentroidAggregator& other) { + centroid_ += other.centroid_; +} + +S2Point CentroidAggregator::Finalize() { + if (centroid_.Norm2() > 0) { + return centroid_.Normalize(); + } else { + return centroid_; + } +} + +void S2ConvexHullAggregator::Add(const S2Geography& geog) { + if (geog.dimension() == 0) { + auto point_ptr = dynamic_cast(&geog); + if (point_ptr != nullptr) { + for (const auto& point : point_ptr->Points()) { + query_.AddPoint(point); + } + } else { + keep_alive_.push_back(s2_rebuild(geog, GlobalOptions())); + Add(*keep_alive_.back()); + } + + return; + } + + if (geog.dimension() == 1) { + auto poly_ptr = dynamic_cast(&geog); + if (poly_ptr != nullptr) { + for (const auto& polyline : poly_ptr->Polylines()) { + query_.AddPolyline(*polyline); + } + } else { + keep_alive_.push_back(s2_rebuild(geog, GlobalOptions())); + Add(*keep_alive_.back()); + } + + return; + } + + if (geog.dimension() == 2) { + auto poly_ptr = dynamic_cast(&geog); + if (poly_ptr != nullptr) { + query_.AddPolygon(*poly_ptr->Polygon()); + } else { + keep_alive_.push_back(s2_rebuild(geog, GlobalOptions())); + Add(*keep_alive_.back()); + } + + return; + } + + auto collection_ptr = dynamic_cast(&geog); + if (collection_ptr != nullptr) { + for (const auto& feature : collection_ptr->Features()) { + Add(*feature); + } + } else { + keep_alive_.push_back(s2_rebuild(geog, GlobalOptions())); + Add(*keep_alive_.back()); + } +} + +std::unique_ptr S2ConvexHullAggregator::Finalize() { + auto polygon = absl::make_unique(); + polygon->Init(query_.GetConvexHull()); + return absl::make_unique(std::move(polygon)); +} + +} diff --git a/src/s2-geography/accessors-geog.hpp b/src/s2-geography/accessors-geog.hpp new file mode 100644 index 00000000..91dc83ed --- /dev/null +++ b/src/s2-geography/accessors-geog.hpp @@ -0,0 +1,35 @@ + +#pragma once + +#include "s2/s2convex_hull_query.h" + +#include "geography.hpp" +#include "aggregator.hpp" + +namespace s2geography { + +S2Point s2_centroid(const S2Geography& geog); +std::unique_ptr s2_boundary(const S2Geography& geog); +std::unique_ptr s2_convex_hull(const S2Geography& geog); + +class CentroidAggregator: public Aggregator { +public: + void Add(const S2Geography& geog); + void Merge(const CentroidAggregator& other); + S2Point Finalize(); + +private: + S2Point centroid_; +}; + +class S2ConvexHullAggregator: public Aggregator> { +public: + void Add(const S2Geography& geog); + std::unique_ptr Finalize(); + +private: + S2ConvexHullQuery query_; + std::vector> keep_alive_; +}; + +} diff --git a/src/s2-geography/accessors.cpp b/src/s2-geography/accessors.cpp new file mode 100644 index 00000000..1ea18d31 --- /dev/null +++ b/src/s2-geography/accessors.cpp @@ -0,0 +1,270 @@ + +#include "geography.hpp" +#include "accessors.hpp" +#include "build.hpp" + +namespace s2geography { + +bool s2_is_collection(const PolygonGeography& geog) { + int num_outer_loops = 0; + for (int i = 0; i < geog.Polygon()->num_loops(); i++) { + S2Loop* loop = geog.Polygon()->loop(i); + num_outer_loops += loop->depth() == 0; + if (num_outer_loops > 1) { + return true; + } + } + + return false; +} + +bool s2_is_collection(const S2Geography& geog) { + int dimension = s2_dimension(geog); + + if (dimension == -1) { + return false; + } + + if (dimension == 0) { + return s2_num_points(geog) > 1; + } + + if (dimension == 1) { + int num_chains = 0; + for (int i = 0; i < geog.num_shapes(); i++) { + std::unique_ptr shape = geog.Shape(i); + num_chains += shape->num_chains(); + if (num_chains > 1) { + return true; + } + } + + return false; + } + + auto polygon_geog_ptr = dynamic_cast(&geog); + if (polygon_geog_ptr != nullptr) { + return s2_is_collection(*polygon_geog_ptr); + } else { + std::unique_ptr built = s2_build_polygon(geog); + return s2_is_collection(*built); + } +} + +int s2_dimension(const S2Geography& geog) { + int dimension = geog.dimension(); + if (dimension != -1) { + return dimension; + } + + for (int i = 0; i < geog.num_shapes(); i++) { + std::unique_ptr shape = geog.Shape(i); + if (shape->dimension() > dimension) { + dimension = shape->dimension(); + } + } + + return dimension; +} + +int s2_num_points(const S2Geography& geog) { + int num_points = 0; + for (int i = 0; i < geog.num_shapes(); i++) { + std::unique_ptr shape = geog.Shape(i); + switch (shape->dimension()) { + case 0: + case 2: + num_points += shape->num_edges(); + break; + case 1: + num_points += shape->num_edges() + shape->num_chains(); + break; + } + } + + return num_points; +} + +bool s2_is_empty(const S2Geography& geog) { + for (int i = 0; i < geog.num_shapes(); i++) { + std::unique_ptr shape = geog.Shape(i); + if (!shape->is_empty()) { + return false; + } + } + + return true; +} + +double s2_area(const PolygonGeography& geog) { + return geog.Polygon()->GetArea(); +} + +double s2_area(const S2GeographyCollection& geog) { + double area = 0; + for (auto& feature: geog.Features()) { + area += s2_area(*feature); + } + return area; +} + +double s2_area(const S2Geography& geog) { + if (s2_dimension(geog) != 2) { + return 0; + } + + auto polygon_geog_ptr = dynamic_cast(&geog); + if (polygon_geog_ptr != nullptr) { + return s2_area(*polygon_geog_ptr); + } + + auto collection_geog_ptr = dynamic_cast(&geog); + if (collection_geog_ptr != nullptr) { + return s2_area(*collection_geog_ptr); + } + + std::unique_ptr built = s2_build_polygon(geog); + return s2_area(*built); +} + +double s2_length(const S2Geography& geog) { + double length = 0; + + if (s2_dimension(geog) == 1) { + for (int i = 0; i < geog.num_shapes(); i++) { + std::unique_ptr shape = geog.Shape(i); + for (int j = 0; j < shape->num_edges(); j++) { + S2Shape::Edge e = shape->edge(j); + S1ChordAngle angle(e.v0, e.v1); + length += angle.radians(); + } + } + } + + return length; +} + +double s2_perimeter(const S2Geography& geog) { + double length = 0; + + if (s2_dimension(geog) == 2) { + for (int i = 0; i < geog.num_shapes(); i++) { + std::unique_ptr shape = geog.Shape(i); + for (int j = 0; j < shape->num_edges(); j++) { + S2Shape::Edge e = shape->edge(j); + S1ChordAngle angle(e.v0, e.v1); + length += angle.radians(); + } + } + } + + return length; +} + +double s2_x(const S2Geography& geog) { + double out = NAN; + for (int i = 0; i < geog.num_shapes(); i++) { + std::unique_ptr shape = geog.Shape(i); + if (shape->dimension() == 0 && shape->num_edges() == 1 && std::isnan(out)) { + S2LatLng pt(shape->edge(0).v0); + out = pt.lng().degrees(); + } else if (shape->dimension() == 0 && shape->num_edges() == 1) { + return NAN; + } + } + + return out; +} + +double s2_y(const S2Geography& geog) { + double out = NAN; + for (int i = 0; i < geog.num_shapes(); i++) { + std::unique_ptr shape = geog.Shape(i); + if (shape->dimension() == 0 && shape->num_edges() == 1 && std::isnan(out)) { + S2LatLng pt(shape->edge(0).v0); + out = pt.lat().degrees(); + } else if (shape->dimension() == 0 && shape->num_edges() == 1) { + return NAN; + } + } + + return out; +} + +bool s2_find_validation_error(const PolylineGeography& geog, S2Error* error) { + for (const auto& polyline: geog.Polylines()) { + if (polyline->FindValidationError(error)) { + return true; + } + } + + return false; +} + +bool s2_find_validation_error(const PolygonGeography& geog, S2Error* error) { + return geog.Polygon()->FindValidationError(error); +} + +bool s2_find_validation_error(const S2GeographyCollection& geog, S2Error* error) { + for (const auto& feature: geog.Features()) { + if (s2_find_validation_error(*feature, error)) { + return true; + } + } + + return false; +} + +bool s2_find_validation_error(const S2Geography& geog, S2Error* error) { + if (geog.dimension() == 0) { + error->Clear(); + return false; + } + + if (geog.dimension() == 1) { + auto poly_ptr = dynamic_cast(&geog); + if (poly_ptr != nullptr) { + return s2_find_validation_error(*poly_ptr, error); + } else { + try { + auto poly = s2_build_polyline(geog); + return s2_find_validation_error(*poly, error); + } catch (Exception& e) { + error->Init(S2Error::INTERNAL, "%s", e.what()); + return true; + } + } + } + + if (geog.dimension() == 2) { + auto poly_ptr = dynamic_cast(&geog); + if (poly_ptr != nullptr) { + return s2_find_validation_error(*poly_ptr, error); + } else { + try { + auto poly = s2_build_polygon(geog); + return s2_find_validation_error(*poly, error); + } catch (Exception& e) { + error->Init(S2Error::INTERNAL, "%s", e.what()); + return true; + } + } + } + + auto collection_ptr = dynamic_cast(&geog); + if (collection_ptr != nullptr) { + return s2_find_validation_error(*collection_ptr, error); + } else { + try { + auto collection = s2_build_polygon(geog); + return s2_find_validation_error(*collection, error); + } catch (Exception& e) { + error->Init(S2Error::INTERNAL, "%s", e.what()); + return true; + } + } + + throw Exception("s2_find_validation() error not implemented for this geography type"); +} + +} diff --git a/src/s2-geography/accessors.hpp b/src/s2-geography/accessors.hpp new file mode 100644 index 00000000..5af41058 --- /dev/null +++ b/src/s2-geography/accessors.hpp @@ -0,0 +1,19 @@ + +#pragma once + +#include "geography.hpp" + +namespace s2geography { + +bool s2_is_collection(const S2Geography& geog); +int s2_dimension(const S2Geography& geog); +int s2_num_points(const S2Geography& geog); +bool s2_is_empty(const S2Geography& geog); +double s2_area(const S2Geography& geog); +double s2_length(const S2Geography& geog); +double s2_perimeter(const S2Geography& geog); +double s2_x(const S2Geography& geog); +double s2_y(const S2Geography& geog); +bool s2_find_validation_error(const S2Geography& geog, S2Error* error); + +} diff --git a/src/s2-geography/aggregator.hpp b/src/s2-geography/aggregator.hpp new file mode 100644 index 00000000..de8e8221 --- /dev/null +++ b/src/s2-geography/aggregator.hpp @@ -0,0 +1,15 @@ + +#pragma once + +#include "geography.hpp" + +namespace s2geography { + +template +class Aggregator { +public: + virtual void Add(const S2Geography& geog, Params... parameters) = 0; + virtual ReturnType Finalize() = 0; +}; + +} diff --git a/src/s2-geography/build.cpp b/src/s2-geography/build.cpp new file mode 100644 index 00000000..8afcc279 --- /dev/null +++ b/src/s2-geography/build.cpp @@ -0,0 +1,409 @@ + +#include "s2/s2boolean_operation.h" +#include "s2/s2builder.h" +#include "s2/s2builderutil_s2polygon_layer.h" +#include "s2/s2builderutil_s2polyline_vector_layer.h" +#include "s2/s2builderutil_s2point_vector_layer.h" +#include "s2/s2builderutil_closed_set_normalizer.h" + +#include "geography.hpp" +#include "build.hpp" +#include "accessors.hpp" + +namespace s2geography { + +std::unique_ptr s2_geography_from_layers(std::vector points, + std::vector> polylines, + std::unique_ptr polygon, + GlobalOptions::OutputAction point_layer_action, + GlobalOptions::OutputAction polyline_layer_action, + GlobalOptions::OutputAction polygon_layer_action) { + // count non-empty dimensions + bool has_polygon = !polygon->is_empty(); + bool has_polylines = polylines.size() > 0; + bool has_points = points.size() > 0; + + // use the requstested dimensions to produce the right kind of EMTPY + bool include_polygon = polygon_layer_action == GlobalOptions::OUTPUT_ACTION_INCLUDE; + bool include_polylines = polyline_layer_action == GlobalOptions::OUTPUT_ACTION_INCLUDE; + bool include_points = point_layer_action == GlobalOptions::OUTPUT_ACTION_INCLUDE; + + if (has_polygon && polygon_layer_action == GlobalOptions::OUTPUT_ACTION_ERROR) { + throw Exception("Output contained unexpected polygon"); + } else if (has_polygon && polygon_layer_action == GlobalOptions::OUTPUT_ACTION_IGNORE) { + has_polygon = false; + } + + if (has_polylines && polyline_layer_action == GlobalOptions::OUTPUT_ACTION_ERROR) { + throw Exception("Output contained unexpected polylines"); + } else if (has_polylines && polyline_layer_action == GlobalOptions::OUTPUT_ACTION_IGNORE) { + has_polylines = false; + } + + if (has_points && point_layer_action == GlobalOptions::OUTPUT_ACTION_ERROR) { + throw Exception("Output contained unexpected points"); + } else if (has_points && point_layer_action == GlobalOptions::OUTPUT_ACTION_IGNORE) { + has_points = false; + } + + int non_empty_dimensions = has_polygon + has_polylines + has_points; + int included_dimensions = include_polygon + include_polylines + include_points; + + // return mixed dimension output + if (non_empty_dimensions > 1) { + std::vector> features; + + if (has_points) { + features.push_back(absl::make_unique(std::move(points))); + } + + if (has_polylines) { + features.push_back(absl::make_unique(std::move(polylines))); + } + + if (has_polygon) { + features.push_back(absl::make_unique(std::move(polygon))); + } + + return absl::make_unique(std::move(features)); + } + + // return single dimension output + if (has_polygon || (included_dimensions == 1 && include_polygon)) { + return absl::make_unique(std::move(polygon)); + } else if (has_polylines || (included_dimensions == 1 && include_polylines)) { + return absl::make_unique(std::move(polylines)); + } else if (has_points || (included_dimensions == 1 && include_points)) { + return absl::make_unique(std::move(points)); + } else { + return absl::make_unique(); + } +} + +std::unique_ptr s2_boolean_operation(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + S2BooleanOperation::OpType op_type, + const GlobalOptions& options) { + + // Create the data structures that will contain the output. + std::vector points; + std::vector> polylines; + std::unique_ptr polygon = absl::make_unique(); + + s2builderutil::LayerVector layers(3); + layers[0] = absl::make_unique(&points, options.point_layer); + layers[1] = absl::make_unique(&polylines, options.polyline_layer); + layers[2] = absl::make_unique(polygon.get(), options.polygon_layer); + + // specify the boolean operation + S2BooleanOperation op( + op_type, + // Normalizing the closed set here is required for line intersections + // to work in the same way as GEOS + s2builderutil::NormalizeClosedSet(std::move(layers)), + options.boolean_operation + ); + + // do the boolean operation, build layers, and check for errors + S2Error error; + if (!op.Build(geog1.ShapeIndex(), geog2.ShapeIndex(), &error)) { + throw Exception(error.text()); + } + + // construct output + return s2_geography_from_layers( + std::move(points), + std::move(polylines), + std::move(polygon), + options.point_layer_action, + options.polyline_layer_action, + options.polygon_layer_action + ); +} + +std::unique_ptr s2_unary_union(const PolygonGeography& geog, + const GlobalOptions& options) { + // A geography with invalid loops won't work with the S2BooleanOperation + // we will use to accumulate (i.e., union) valid polygons, + // so we need to rebuild each loop as its own polygon, + // splitting crossed edges along the way. + + // Not exposing these options as an argument (except snap function) + // because a particular combiation of them is required for this to work + S2Builder::Options builder_options; + builder_options.set_split_crossing_edges(true); + builder_options.set_snap_function(options.boolean_operation.snap_function()); + s2builderutil::S2PolygonLayer::Options layer_options; + layer_options.set_edge_type(S2Builder::EdgeType::UNDIRECTED); + layer_options.set_validate(false); + + // Rebuild all loops as polygons using the S2Builder() + std::vector> loops; + for (int i = 0; i < geog.Polygon()->num_loops(); i++) { + std::unique_ptr loop = absl::make_unique(); + S2Builder builder(builder_options); + builder.StartLayer( + absl::make_unique(loop.get(), layer_options)); + builder.AddShape(S2Loop::Shape(geog.Polygon()->loop(i))); + S2Error error; + if (!builder.Build(&error)) { + throw Exception(error.text()); + } + + // Check if the builder created a polygon whose boundary contained more than + // half the earth (and invert it if so) + if (loop->GetArea() > (2 * M_PI)) { + loop->Invert(); + } + + loops.push_back(std::move(loop)); + } + + // Accumulate the union of outer loops (but difference of inner loops) + std::unique_ptr accumulated_polygon = absl::make_unique(); + for (int i = 0; i < geog.Polygon()->num_loops(); i++) { + std::unique_ptr polygon_result = absl::make_unique(); + + // Use original nesting to suggest if this loop should be unioned or diffed. + // For valid polygons loops are arranged such that the biggest loop is on the outside + // followed by holes such that the below strategy should work (since we are + // just iterating along the original loop structure) + if ((geog.Polygon()->loop(i)->depth() % 2) == 0) { + polygon_result->InitToUnion(accumulated_polygon.get(), loops[i].get()); + } else { + polygon_result->InitToDifference(accumulated_polygon.get(), loops[i].get()); + } + + accumulated_polygon.swap(polygon_result); + } + + return absl::make_unique(std::move(accumulated_polygon)); +} + +std::unique_ptr s2_unary_union(const ShapeIndexGeography& geog, + const GlobalOptions& options) { + // complex union only needed when a polygon is involved + bool simple_union_ok = s2_is_empty(geog) || s2_dimension(geog) < 2; + + // valid polygons that are not part of a collection can also use a + // simple union (common) + if (geog.dimension() == 2) { + S2Error validation_error; + if (!s2_find_validation_error(geog, &validation_error)) { + simple_union_ok = true; + } + } + + if (simple_union_ok) { + ShapeIndexGeography empty; + return s2_boolean_operation(geog, empty, S2BooleanOperation::OpType::UNION, options); + } + + if (geog.dimension() == 2) { + // If we've made it here we have an invalid polygon on our hands. + auto poly_ptr = dynamic_cast(&geog); + if (poly_ptr != nullptr) { + return s2_unary_union(*poly_ptr, options); + } else { + auto poly = s2_build_polygon(geog); + return s2_unary_union(*poly, options); + } + } + + throw Exception( + "s2_unary_union() for multidimensional collections not implemented"); +} + +std::unique_ptr s2_rebuild(const S2Geography& geog, + const GlobalOptions& options, + GlobalOptions::OutputAction point_layer_action, + GlobalOptions::OutputAction polyline_layer_action, + GlobalOptions::OutputAction polygon_layer_action) { + // create the builder + S2Builder builder(options.builder); + + // create the data structures that will contain the output + std::vector points; + std::vector> polylines; + std::unique_ptr polygon = absl::make_unique(); + + // add shapes to the layer with the appropriate dimension + builder.StartLayer( + absl::make_unique(&points, options.point_layer) + ); + for (int i = 0; i < geog.num_shapes(); i++) { + auto shape = geog.Shape(i); + if (shape->dimension() == 0) { + builder.AddShape(*shape); + } + } + + builder.StartLayer( + absl::make_unique(&polylines, options.polyline_layer) + ); + for (int i = 0; i < geog.num_shapes(); i++) { + auto shape = geog.Shape(i); + if (shape->dimension() == 1) { + builder.AddShape(*shape); + } + } + + builder.StartLayer( + absl::make_unique(polygon.get(), options.polygon_layer) + ); + for (int i = 0; i < geog.num_shapes(); i++) { + auto shape = geog.Shape(i); + if (shape->dimension() == 2) { + builder.AddShape(*shape); + } + } + + // build the output + S2Error error; + if (!builder.Build(&error)) { + throw Exception(error.text()); + } + + // construct output + return s2_geography_from_layers( + std::move(points), + std::move(polylines), + std::move(polygon), + options.point_layer_action, + options.polyline_layer_action, + options.polygon_layer_action + ); +} + +std::unique_ptr s2_rebuild(const S2Geography& geog, + const GlobalOptions& options) { + return s2_rebuild( + geog, + options, + options.point_layer_action, + options.polyline_layer_action, + options.polygon_layer_action); +} + +std::unique_ptr s2_build_point(const S2Geography& geog) { + std::unique_ptr geog_out = s2_rebuild( + geog, + GlobalOptions(), + GlobalOptions::OutputAction::OUTPUT_ACTION_INCLUDE, + GlobalOptions::OutputAction::OUTPUT_ACTION_ERROR, + GlobalOptions::OutputAction::OUTPUT_ACTION_ERROR); + + return std::unique_ptr( + dynamic_cast(geog_out.release())); +} + + +std::unique_ptr s2_build_polyline(const S2Geography& geog) { + std::unique_ptr geog_out = s2_rebuild( + geog, + GlobalOptions(), + GlobalOptions::OutputAction::OUTPUT_ACTION_ERROR, + GlobalOptions::OutputAction::OUTPUT_ACTION_INCLUDE, + GlobalOptions::OutputAction::OUTPUT_ACTION_ERROR); + + return std::unique_ptr( + dynamic_cast(geog_out.release())); +} + + +std::unique_ptr s2_build_polygon(const S2Geography& geog) { + std::unique_ptr geog_out = s2_rebuild( + geog, + GlobalOptions(), + GlobalOptions::OutputAction::OUTPUT_ACTION_ERROR, + GlobalOptions::OutputAction::OUTPUT_ACTION_ERROR, + GlobalOptions::OutputAction::OUTPUT_ACTION_INCLUDE); + + return std::unique_ptr( + dynamic_cast(geog_out.release())); +} + + +void RebuildAggregator::Add(const S2Geography& geog) { + index_.Add(geog); +} + +std::unique_ptr RebuildAggregator::Finalize() { + return s2_rebuild(index_, options_); +} + +void S2CoverageUnionAggregator::Add(const S2Geography& geog) { + index_.Add(geog); +} + +std::unique_ptr S2CoverageUnionAggregator::Finalize() { + ShapeIndexGeography empty_index_; + return s2_boolean_operation(index_, empty_index_, S2BooleanOperation::OpType::UNION, options_); +} + +void S2UnionAggregator::Add(const S2Geography& geog) { + if (geog.dimension() == 0 || geog.dimension() == 1) { + root_.index1.Add(geog); + return; + } + + if (other_.size() == 0) { + other_.push_back(absl::make_unique()); + other_.back()->index1.Add(geog); + return; + } + + Node* last = other_.back().get(); + if (last->index1.num_shapes() == 0) { + last->index1.Add(geog); + } else if (last->index2.num_shapes() == 0) { + last->index2.Add(geog); + } else { + other_.push_back(absl::make_unique()); + other_.back()->index1.Add(geog); + } +} + +std::unique_ptr S2UnionAggregator::Node::Merge(const GlobalOptions& options) { + return s2_boolean_operation( + index1, + index2, + S2BooleanOperation::OpType::UNION, + options); +} + +std::unique_ptr S2UnionAggregator::Finalize() { + for (int j = 0; j < 100; j++) { + if (other_.size() <= 1) { + break; + } + + for (int64_t i = static_cast(other_.size()) - 1; i >= 1; i = i - 2) { + // merge other_[i] with other_[i - 1] + std::unique_ptr merged = other_[i]->Merge(options_); + std::unique_ptr merged_prev = other_[i - 1]->Merge(options_); + + // erase the last two nodes + other_.erase(other_.begin() + i - 1, other_.begin() + i + 1); + + // ..and replace it with a single node + other_.push_back(absl::make_unique()); + other_.back()->index1.Add(*merged); + other_.back()->index2.Add(*merged_prev); + + // making sure to keep the underlying data alive + other_.back()->data.push_back(std::move(merged)); + other_.back()->data.push_back(std::move(merged_prev)); + } + } + + if (other_.size() == 0) { + return root_.Merge(options_); + } else { + std::unique_ptr merged = other_[0]->Merge(options_); + root_.index2.Add(*merged); + return root_.Merge(options_); + } +} + +} diff --git a/src/s2-geography/build.hpp b/src/s2-geography/build.hpp new file mode 100644 index 00000000..fde82254 --- /dev/null +++ b/src/s2-geography/build.hpp @@ -0,0 +1,97 @@ + +#pragma once + +#include "s2/s2builderutil_s2polygon_layer.h" +#include "s2/s2builderutil_s2polyline_vector_layer.h" +#include "s2/s2builderutil_s2point_vector_layer.h" + +#include "geography.hpp" +#include "aggregator.hpp" + +namespace s2geography { + +class GlobalOptions { +public: + enum OutputAction { + OUTPUT_ACTION_INCLUDE, + OUTPUT_ACTION_IGNORE, + OUTPUT_ACTION_ERROR + }; + + GlobalOptions() + : point_layer_action(OUTPUT_ACTION_INCLUDE), + polyline_layer_action(OUTPUT_ACTION_INCLUDE), + polygon_layer_action(OUTPUT_ACTION_INCLUDE) {} + + S2BooleanOperation::Options boolean_operation; + S2Builder::Options builder; + + s2builderutil::S2PointVectorLayer::Options point_layer; + s2builderutil::S2PolylineVectorLayer::Options polyline_layer; + s2builderutil::S2PolygonLayer::Options polygon_layer; + OutputAction point_layer_action; + OutputAction polyline_layer_action; + OutputAction polygon_layer_action; +}; + +std::unique_ptr s2_boolean_operation(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + S2BooleanOperation::OpType op_type, + const GlobalOptions& options); + +std::unique_ptr s2_unary_union(const ShapeIndexGeography& geog, + const GlobalOptions& options); + +std::unique_ptr s2_rebuild(const S2Geography& geog, + const GlobalOptions& options); + +std::unique_ptr s2_build_point(const S2Geography& geog); + +std::unique_ptr s2_build_polyline(const S2Geography& geog); + +std::unique_ptr s2_build_polygon(const S2Geography& geog); + +class RebuildAggregator: public Aggregator> { +public: + RebuildAggregator(const GlobalOptions& options): options_(options) {} + void Add(const S2Geography& geog); + std::unique_ptr Finalize(); + +private: + GlobalOptions options_; + ShapeIndexGeography index_; +}; + +class S2CoverageUnionAggregator: public Aggregator> { +public: + S2CoverageUnionAggregator(const GlobalOptions& options): options_(options) {} + + void Add(const S2Geography& geog); + std::unique_ptr Finalize(); + +private: + GlobalOptions options_; + ShapeIndexGeography index_; +}; + +class S2UnionAggregator: public Aggregator> { +public: + S2UnionAggregator(const GlobalOptions& options): options_(options) {} + void Add(const S2Geography& geog); + std::unique_ptr Finalize(); + +private: + class Node { + public: + ShapeIndexGeography index1; + ShapeIndexGeography index2; + std::vector> data; + std::unique_ptr Merge(const GlobalOptions& options); + }; + + GlobalOptions options_; + Node root_; + std::vector> other_; +}; + +} diff --git a/src/s2-geography/constructor.hpp b/src/s2-geography/constructor.hpp new file mode 100644 index 00000000..491d7712 --- /dev/null +++ b/src/s2-geography/constructor.hpp @@ -0,0 +1,346 @@ + +#pragma once + +#include + +#include "geography.hpp" + +namespace s2geography { + +namespace util { + +enum GeometryType { + GEOMETRY_TYPE_UNKNOWN = 0, + POINT = 1, + LINESTRING = 2, + POLYGON = 3, + MULTIPOINT = 4, + MULTILINESTRING = 5, + MULTIPOLYGON = 6, + GEOMETRYCOLLECTION = 7 +}; + +} + +class Constructor { +public: + + class Options { + public: + Options() : oriented_(false), check_(true) {} + bool oriented() { return oriented_; } + void set_oriented(bool oriented) { oriented_ = oriented; } + bool check() { return check_; } + void set_check(bool check) { check_ = check; } + + private: + bool oriented_; + bool check_; + }; + + Constructor(const Options& options): options_(options) {} + + virtual ~Constructor() {} + + Options* mutable_options() { return &options_; } + + virtual void geom_start(util::GeometryType geometry_type, int64_t size) {} + virtual void ring_start(int32_t size) {} + + virtual void coords(const double* coord, int64_t n, int32_t coord_size) { + for (int64_t i = 0; i < n; i++) { + S2LatLng pt = S2LatLng::FromDegrees(coord[i * coord_size + 1], coord[i * coord_size]); + points_.push_back(pt.Normalized().ToPoint()); + } + } + + virtual void ring_end() {} + virtual void geom_end() {} + virtual std::unique_ptr finish() = 0; + +protected: + std::vector points_; + Options options_; +}; + +class PointConstructor: public Constructor { +public: + + PointConstructor(): Constructor(Options()) {} + + void geom_start(util::GeometryType geometry_type, int64_t size) { + if (size != 0 && + geometry_type != util::GeometryType::POINT && + geometry_type != util::GeometryType::MULTIPOINT && + geometry_type != util::GeometryType::GEOMETRYCOLLECTION) { + throw Exception( + "PointConstructor input must be empty, point, multipoint, or collection"); + } + + if (size > 0) { + points_.reserve(points_.size() + size); + } + } + + void coords(const double* coord, int64_t n, int32_t coord_size) { + for (int64_t i = 0; i < n; i++) { + if (coord_empty(coord + (i * coord_size), coord_size)) { + continue; + } + + S2LatLng pt = S2LatLng::FromDegrees(coord[i * coord_size + 1], coord[i * coord_size]); + points_.push_back(pt.ToPoint()); + } + } + + std::unique_ptr finish() { + auto result = absl::make_unique(std::move(points_)); + points_.clear(); + return std::unique_ptr(result.release()); + } + +private: + bool coord_empty(const double* coord, int32_t coord_size) { + for (int32_t i = 0; i < coord_size; i++) { + if (!std::isnan(coord[i])) { + return false; + } + } + + return true; + } +}; + +class PolylineConstructor: public Constructor { +public: + PolylineConstructor(const Options& options): Constructor(options) {} + + void geom_start(util::GeometryType geometry_type, int64_t size) { + if (size != 0 && + geometry_type != util::GeometryType::LINESTRING && + geometry_type != util::GeometryType::MULTILINESTRING && + geometry_type != util::GeometryType::GEOMETRYCOLLECTION) { + throw Exception( + "PolylineConstructor input must be empty, linestring, multilinestring, or collection"); + } + + if (size > 0 && geometry_type == util::GeometryType::LINESTRING) { + points_.reserve(size); + } + } + + void geom_end() { + if (points_.size() > 0) { + auto polyline = absl::make_unique(); + polyline->Init(std::move(points_)); + + if (options_.check() && !polyline->IsValid()) { + polyline->FindValidationError(&error_); + throw Exception(error_.text()); + } + + polylines_.push_back(std::move(polyline)); + points_.clear(); + } + } + + std::unique_ptr finish() { + std::unique_ptr result; + + if (polylines_.size() > 0) { + result = absl::make_unique(std::move(polylines_)); + polylines_.clear(); + } else { + result = absl::make_unique(); + } + + return std::unique_ptr(result.release()); + } + +private: + std::vector> polylines_; + S2Error error_; +}; + +class PolygonConstructor: public Constructor { +public: + PolygonConstructor(const Options& options): Constructor(options) {} + + void ring_start(int32_t size) { + points_.clear(); + if (size > 0) { + points_.reserve(size); + } + } + + void ring_end() { + if (points_.size() == 0) { + return; + } + + // S2Loop is open instead of closed + points_.pop_back(); + auto loop = absl::make_unique(); + loop->set_s2debug_override(S2Debug::DISABLE); + loop->Init(std::move(points_)); + + if (!options_.oriented()) { + loop->Normalize(); + } + + if (options_.check() && !loop->IsValid()) { + std::stringstream err; + err << "Loop " << (loops_.size()) << " is not valid: "; + loop->FindValidationError(&error_); + err << error_.text(); + throw Exception(err.str()); + } + + loops_.push_back(std::move(loop)); + points_.clear(); + } + + std::unique_ptr finish() { + auto polygon = absl::make_unique(); + polygon->set_s2debug_override(S2Debug::DISABLE); + if (options_.oriented()) { + polygon->InitOriented(std::move(loops_)); + } else { + polygon->InitNested(std::move(loops_)); + } + + loops_.clear(); + + if (options_.check() && !polygon->IsValid()) { + polygon->FindValidationError(&error_); + throw Exception(error_.text()); + } + + auto result = absl::make_unique(std::move(polygon)); + return std::unique_ptr(result.release()); + } + +private: + std::vector> loops_; + S2Error error_; +}; + +class CollectionConstructor: public Constructor { +public: + CollectionConstructor(const Options& options): + Constructor(options), + polyline_constructor_(options), + polygon_constructor_(options), + collection_constructor_(nullptr), + level_(0) {} + + void geom_start(util::GeometryType geometry_type, int64_t size) { + level_++; + if (level_ == 1 && geometry_type == util::GeometryType::GEOMETRYCOLLECTION) { + active_constructor_ = nullptr; + return; + } + + if (active_constructor_ != nullptr) { + active_constructor_->geom_start(geometry_type, size); + return; + } + + switch (geometry_type) { + case util::GeometryType::POINT: + case util::GeometryType::MULTIPOINT: + active_constructor_ = &point_constructor_; + break; + case util::GeometryType::LINESTRING: + case util::GeometryType::MULTILINESTRING: + active_constructor_ = &polyline_constructor_; + break; + case util::GeometryType::POLYGON: + case util::GeometryType::MULTIPOLYGON: + active_constructor_ = &polygon_constructor_; + break; + case util::GeometryType::GEOMETRYCOLLECTION: + this->collection_constructor_ = absl::make_unique(options_); + this->active_constructor_ = this->collection_constructor_.get(); + break; + default: + throw Exception("CollectionConstructor: unsupported geometry type"); + } + + active_constructor_->geom_start(geometry_type, size); + } + + void ring_start(int32_t size) { + active_constructor_->ring_start(size); + } + + void coords(const double* coord, int64_t n, int32_t coord_size) { + active_constructor_->coords(coord, n, coord_size); + } + + void ring_end() { + active_constructor_->ring_end(); + } + + void geom_end() { + level_--; + + if (level_ >= 1) { + active_constructor_->geom_end(); + } + + if (level_ == 1) { + auto feature = active_constructor_->finish(); + features_.push_back(std::move(feature)); + active_constructor_ = nullptr; + } + } + + std::unique_ptr finish() { + auto result = absl::make_unique(std::move(features_)); + features_.clear(); + return std::unique_ptr(result.release()); + } + +private: + PointConstructor point_constructor_; + PolylineConstructor polyline_constructor_; + PolygonConstructor polygon_constructor_; + std::unique_ptr collection_constructor_; + + +protected: + Constructor* active_constructor_; + int level_; + std::vector> features_; +}; + +class VectorConstructor: public CollectionConstructor { +public: + VectorConstructor(const Options& options): CollectionConstructor(options) {} + + void start_feature() { + active_constructor_ = nullptr; + level_ = 0; + features_.clear(); + geom_start(util::GeometryType::GEOMETRYCOLLECTION, 1); + } + + std::unique_ptr finish_feature() { + geom_end(); + + if (features_.empty()) { + return absl::make_unique(); + } else { + std::unique_ptr feature = std::move(features_.back()); + if (feature.get() == nullptr) { + throw Exception("finish_feature() generated nullptr"); + } + + features_.pop_back(); + return feature; + } + } +}; + +} diff --git a/src/s2-geography/coverings.cpp b/src/s2-geography/coverings.cpp new file mode 100644 index 00000000..eb8dd8d9 --- /dev/null +++ b/src/s2-geography/coverings.cpp @@ -0,0 +1,78 @@ + +#include "s2/s2region_coverer.h" +#include "s2/s2shape_index_buffered_region.h" + +#include "geography.hpp" +#include "coverings.hpp" +#include "accessors.hpp" +#include "accessors-geog.hpp" + +namespace s2geography { + +S2Point s2_point_on_surface(const S2Geography& geog, S2RegionCoverer& coverer) { + if (s2_is_empty(geog)) { + return S2Point(); + } + + int dimension = s2_dimension(geog); + if (dimension == 2) { + std::unique_ptr region = geog.Region(); + S2CellUnion covering = coverer.GetInteriorCovering(*region); + + // Take center of cell with smallest level (biggest) + int min_level = 31; + S2Point pt; + for (const S2CellId& id : covering) { + if (id.level() < min_level) { + // Already normalized + pt = id.ToPoint(); + min_level = id.level(); + } + } + + return pt; + } + + if (dimension == 0) { + // For point, return point closest to centroid + S2Point centroid = s2_centroid(geog); + + S1Angle nearest_dist = S1Angle::Infinity(); + S1Angle dist; + S2Point closest_pt; + for (int i = 0; i < geog.num_shapes(); i++) { + auto shape = geog.Shape(i); + for (int j = 0; j < shape->num_edges(); j++) { + S2Shape::Edge e = shape->edge(j); + dist = S1Angle(e.v0, centroid); + if (dist < nearest_dist) { + nearest_dist = dist; + closest_pt = e.v0; + } + } + } + + return closest_pt; + } + + throw Exception("s2_point_on_surface() not implemented for polyline"); +} + +void s2_covering(const S2Geography& geog, std::vector* covering, + S2RegionCoverer& coverer) { + coverer.GetCovering(*geog.Region(), covering); +} + +void s2_interior_covering(const S2Geography& geog, std::vector* covering, + S2RegionCoverer& coverer) { + coverer.GetInteriorCovering(*geog.Region(), covering); +} + +void s2_covering_buffered(const ShapeIndexGeography& geog, double distance_radians, + std::vector* covering, + S2RegionCoverer& coverer) { + S2ShapeIndexBufferedRegion region(&geog.ShapeIndex(), S1ChordAngle::Radians(distance_radians)); + coverer.GetCovering(region, covering); +} + +} diff --git a/src/s2-geography/coverings.hpp b/src/s2-geography/coverings.hpp new file mode 100644 index 00000000..d9d6ceb3 --- /dev/null +++ b/src/s2-geography/coverings.hpp @@ -0,0 +1,19 @@ + +#pragma once + +#include "s2/s2region_coverer.h" + +#include "geography.hpp" + +namespace s2geography { + +S2Point s2_point_on_surface(const S2Geography& geog, S2RegionCoverer& coverer); +void s2_covering(const S2Geography& geog, std::vector* covering, + S2RegionCoverer& coverer); +void s2_interior_covering(const S2Geography& geog, std::vector* covering, + S2RegionCoverer& coverer); +void s2_covering_buffered(const ShapeIndexGeography& geog, double distance_radians, + std::vector* covering, + S2RegionCoverer& coverer); + +} diff --git a/src/s2-geography/distance.cpp b/src/s2-geography/distance.cpp new file mode 100644 index 00000000..24fd4509 --- /dev/null +++ b/src/s2-geography/distance.cpp @@ -0,0 +1,69 @@ + +#include "s2/s2closest_edge_query.h" +#include "s2/s2furthest_edge_query.h" + +#include "geography.hpp" +#include "distance.hpp" + +namespace s2geography { + +double s2_distance(const ShapeIndexGeography& geog1, const ShapeIndexGeography& geog2) { + S2ClosestEdgeQuery query(&geog1.ShapeIndex()); + S2ClosestEdgeQuery::ShapeIndexTarget target(&geog2.ShapeIndex()); + + const auto& result = query.FindClosestEdge(&target); + + S1ChordAngle angle = result.distance(); + return angle.ToAngle().radians(); +} + +double s2_max_distance(const ShapeIndexGeography& geog1, const ShapeIndexGeography& geog2) { + S2FurthestEdgeQuery query(&geog1.ShapeIndex()); + S2FurthestEdgeQuery::ShapeIndexTarget target(&geog2.ShapeIndex()); + + const auto& result = query.FindFurthestEdge(&target); + + S1ChordAngle angle = result.distance(); + return angle.ToAngle().radians(); +} + +S2Point s2_closest_point(const ShapeIndexGeography& geog1, const ShapeIndexGeography& geog2) { + return s2_minimum_clearance_line_between(geog1, geog2).first; +} + +std::pair s2_minimum_clearance_line_between( + const ShapeIndexGeography& geog1, const ShapeIndexGeography& geog2) { + S2ClosestEdgeQuery query1(&geog1.ShapeIndex()); + query1.mutable_options()->set_include_interiors(false); + S2ClosestEdgeQuery::ShapeIndexTarget target(&geog2.ShapeIndex()); + + const auto& result1 = query1.FindClosestEdge(&target); + + if (result1.edge_id() == -1) { + return std::pair(S2Point(0, 0, 0), S2Point(0, 0, 0)); + } + + // Get the edge from index1 (edge1) that is closest to index2. + S2Shape::Edge edge1 = query1.GetEdge(result1); + + // Now find the edge from index2 (edge2) that is closest to edge1. + S2ClosestEdgeQuery query2(&geog2.ShapeIndex()); + query2.mutable_options()->set_include_interiors(false); + S2ClosestEdgeQuery::EdgeTarget target2(edge1.v0, edge1.v1); + auto result2 = query2.FindClosestEdge(&target2); + + // what if result2 has no edges? + if (result2.is_interior()) { + throw Exception("S2ClosestEdgeQuery result is interior!"); + } + + S2Shape::Edge edge2 = query2.GetEdge(result2); + + // Find the closest point pair on edge1 and edge2. + return S2::GetEdgePairClosestPoints( + edge1.v0, edge1.v1, + edge2.v0, edge2.v1 + ); +} + +} diff --git a/src/s2-geography/distance.hpp b/src/s2-geography/distance.hpp new file mode 100644 index 00000000..cd257298 --- /dev/null +++ b/src/s2-geography/distance.hpp @@ -0,0 +1,14 @@ + +#pragma once + +#include "geography.hpp" + +namespace s2geography { + +double s2_distance(const ShapeIndexGeography& geog1, const ShapeIndexGeography& geog2); +double s2_max_distance(const ShapeIndexGeography& geog1, const ShapeIndexGeography& geog2); +S2Point s2_closest_point(const ShapeIndexGeography& geog1, const ShapeIndexGeography& geog2); +std::pair s2_minimum_clearance_line_between( + const ShapeIndexGeography& geog1, const ShapeIndexGeography& geog2); + +} diff --git a/src/s2-geography/geography.cpp b/src/s2-geography/geography.cpp new file mode 100644 index 00000000..04aa4e25 --- /dev/null +++ b/src/s2-geography/geography.cpp @@ -0,0 +1,172 @@ + +#include "s2/s2shape.h" +#include "s2/s2region.h" +#include "s2/s2region_union.h" + +#include "s2/s2point_region.h" +#include "s2/s2point_vector_shape.h" +#include "s2/s2polyline.h" +#include "s2/s2polygon.h" + +#include "s2/mutable_s2shape_index.h" +#include "s2/s2shape_index_region.h" +#include "s2/s2shapeutil_coding.h" + +#include "geography.hpp" +using namespace s2geography; + + +// This class is a shim to allow a class to return a std::unique_ptr(), +// which is required by MutableS2ShapeIndex::Add(), without copying the underlying +// data. S2Shape instances do not typically own their data (e.g., S2Polygon::Shape), +// so this does not change the general relationship (that anything returned by +// S2Geography::Shape() is only valid within the scope of the S2Geography). +// Note that this class is also available (but not exposed) in s2/s2shapeutil_coding.cc. +class S2ShapeWrapper: public S2Shape { +public: + S2ShapeWrapper(S2Shape* shape): shape_(shape) {} + int num_edges() const { return shape_->num_edges(); } + Edge edge(int edge_id) const { return shape_->edge(edge_id); } + int dimension() const { return shape_->dimension(); } + ReferencePoint GetReferencePoint() const { return shape_->GetReferencePoint(); } + int num_chains() const { return shape_->num_chains(); } + Chain chain(int chain_id) const { return shape_->chain(chain_id); } + Edge chain_edge(int chain_id, int offset) const { return shape_->chain_edge(chain_id, offset); } + ChainPosition chain_position(int edge_id) const { return shape_->chain_position(edge_id); } + +private: + S2Shape* shape_; +}; + +// Just like the S2ShapeWrapper, the S2RegionWrapper helps reconcile the differences +// in lifecycle expectation between S2 and S2Geography. We often need access to a +// S2Region to generalize algorithms; however, there are some operations that need +// ownership of the region (e.g., the S2RegionUnion). In S2Geography the assumption +// is that anything returned by a S2Geography is only valid for the lifetime of the +// underlying S2Geography. A different design of the algorithms implemented here might +// well make this unnecessary. +class S2RegionWrapper: public S2Region { +public: + S2RegionWrapper(S2Region* region): region_(region) {} + S2Region* Clone() const { return region_->Clone(); } + S2Cap GetCapBound() const { return region_->GetCapBound(); } + S2LatLngRect GetRectBound() const { return region_->GetRectBound(); } + void GetCellUnionBound(std::vector *cell_ids) const { + return region_->GetCellUnionBound(cell_ids); + } + bool Contains(const S2Cell& cell) const { return region_->Contains(cell); } + bool MayIntersect(const S2Cell& cell) const { return region_->MayIntersect(cell); } + bool Contains(const S2Point& p) const { return region_->Contains(p); } + +private: + S2Region* region_; +}; + + +void S2Geography::GetCellUnionBound(std::vector* cell_ids) const { + MutableS2ShapeIndex index; + for (int i = 0; i < num_shapes(); i++) { + index.Add(Shape(i)); + } + + MakeS2ShapeIndexRegion(&index).GetCellUnionBound(cell_ids); +} + +std::unique_ptr PointGeography::Shape(int id) const { + return absl::make_unique(points_); +} + +std::unique_ptr PointGeography::Region() const { + auto region = absl::make_unique(); + for (const S2Point& point: points_) { + region->Add(absl::make_unique(point)); + } + + // because Rtools for R 3.6 on Windows complains about a direct + // return region + return std::unique_ptr(region.release()); +} + +void PointGeography::GetCellUnionBound(std::vector* cell_ids) const { + if (points_.size() < 10) { + for (const S2Point& point: points_) { + cell_ids->push_back(S2CellId(point)); + } + } else { + S2Geography::GetCellUnionBound(cell_ids); + } +} + + +int PolylineGeography::num_shapes() const { return polylines_.size(); } + +std::unique_ptr PolylineGeography::Shape(int id) const { + return absl::make_unique(polylines_[id].get()); +} + +std::unique_ptr PolylineGeography::Region() const { + auto region = absl::make_unique(); + for (const auto& polyline: polylines_) { + region->Add(absl::make_unique(polyline.get())); + } + // because Rtools for R 3.6 on Windows complains about a direct + // return region + return std::unique_ptr(region.release()); +} + +void PolylineGeography::GetCellUnionBound(std::vector* cell_ids) const { + for (const auto& polyline: polylines_) { + polyline->GetCellUnionBound(cell_ids); + } +} + +std::unique_ptr PolygonGeography::Shape(int id) const { + return absl::make_unique(polygon_.get()); +} + +std::unique_ptr PolygonGeography::Region() const { + return absl::make_unique(polygon_.get()); +} + +void PolygonGeography::GetCellUnionBound(std::vector* cell_ids) const { + polygon_->GetCellUnionBound(cell_ids); +} + + +int S2GeographyCollection::num_shapes() const { return total_shapes_; } + +std::unique_ptr S2GeographyCollection::Shape(int id) const { + int sum_shapes_ = 0; + for (int i = 0; i < features_.size(); i++) { + sum_shapes_ += num_shapes_[i]; + if (id < sum_shapes_) { + return features_[i]->Shape(id - sum_shapes_ + num_shapes_[i]); + } + } + + throw Exception("shape id out of bounds"); +} + +std::unique_ptr S2GeographyCollection::Region() const { + auto region = absl::make_unique(); + for (const auto& feature: features_) { + region->Add(feature->Region()); + } + // because Rtools for R 3.6 on Windows complains about a direct + // return region + return std::unique_ptr(region.release()); +} + +int ShapeIndexGeography::num_shapes() const { return shape_index_.num_shape_ids(); } + +std::unique_ptr ShapeIndexGeography::Shape(int id) const { + S2Shape* shape = shape_index_.shape(id); + return std::unique_ptr(new S2ShapeWrapper(shape)); +} + +std::unique_ptr ShapeIndexGeography::Region() const { + auto region = absl::make_unique>(&shape_index_); + // because Rtools for R 3.6 on Windows complains about a direct + // return region + return std::unique_ptr(region.release()); +} diff --git a/src/s2-geography/geography.hpp b/src/s2-geography/geography.hpp new file mode 100644 index 00000000..db3e040a --- /dev/null +++ b/src/s2-geography/geography.hpp @@ -0,0 +1,215 @@ + +#pragma once + +#include +#include +#include + +#include "s2/s2polygon.h" +#include "s2/s2polyline.h" + +namespace s2geography { + +class Exception: public std::runtime_error { +public: + Exception(std::string what): std::runtime_error(what.c_str()) {} +}; + +// An S2Geography is an abstraction of S2 types that is designed to closely match +// the scope of a GEOS Geometry. Its methods are limited to those needed to +// implement C API functions. From an S2 perspective, an S2Geography is +// an S2Region that can be represented by zero or more S2Shape objects. +// Current implementations of S2Geography own their data (i.e., the +// coordinate vectors and underlying S2 objects), however, the interface is +// designed to allow future abstractions where this is not the case. +class S2Geography { +public: + + virtual ~S2Geography() {} + + // Returns 0, 1, or 2 if all Shape()s that are returned will have + // the same dimension (i.e., they are all points, all lines, or + // all polygons). + virtual int dimension() const { + if (num_shapes() == 0) { + return -1; + } + + int dim = Shape(0)->dimension(); + for (int i = 1; i < num_shapes(); i++) { + if (dim != Shape(i)->dimension()) { + return -1; + } + } + + return dim; + } + + // The number of S2Shape objects needed to represent this S2Geography + virtual int num_shapes() const = 0; + + // Returns the given S2Shape (where 0 <= id < num_shapes()). The + // caller retains ownership of the S2Shape but the data pointed to + // by the object requires that the underlying S2Geography outlives + // the returned object. + virtual std::unique_ptr Shape(int id) const = 0; + + // Returns an S2Region that represents the object. The caller retains + // ownership of the S2Region but the data pointed to by the object + // requires that the underlying S2Geography outlives the returned + // object. + virtual std::unique_ptr Region() const = 0; + + // Adds an unnormalized set of S2CellIDs to `cell_ids`. This is intended + // to be faster than using Region().GetCovering() directly and to + // return a small number of cells that can be used to compute a possible + // intersection quickly. + virtual void GetCellUnionBound(std::vector* cell_ids) const; +}; + + +// An S2Geography representing zero or more points using a std::vector +// as the underlying representation. +class PointGeography: public S2Geography { +public: + PointGeography() {} + PointGeography(S2Point point) { points_.push_back(point); } + PointGeography(std::vector points): points_(std::move(points)) {} + + int dimension() const { return 0; } + int num_shapes() const { return 1; } + std::unique_ptr Shape(int id) const; + std::unique_ptr Region() const; + void GetCellUnionBound(std::vector* cell_ids) const; + + const std::vector& Points() const { + return points_; + } + +private: + std::vector points_; +}; + + +// An S2Geography representing zero or more polylines using the S2Polyline class +// as the underlying representation. +class PolylineGeography: public S2Geography { +public: + PolylineGeography() {} + PolylineGeography(std::unique_ptr polyline) { + polylines_.push_back(std::move(polyline)); + } + PolylineGeography(std::vector> polylines): + polylines_(std::move(polylines)) {} + + int dimension() const { return 1; } + int num_shapes() const; + std::unique_ptr Shape(int id) const; + std::unique_ptr Region() const; + void GetCellUnionBound(std::vector* cell_ids) const; + + const std::vector>& Polylines() const { + return polylines_; + } + +private: + std::vector> polylines_; +}; + + +// An S2Geography representing zero or more polygons using the S2Polygon class +// as the underlying representation. Note that a single S2Polygon (from the S2 +// perspective) can represent zero or more polygons (from the simple features +// perspective). +class PolygonGeography: public S2Geography { +public: + PolygonGeography() {} + PolygonGeography(std::unique_ptr polygon): + polygon_(std::move(polygon)) {} + + int dimension() const { return 2; } + int num_shapes() const { return 1; } + std::unique_ptr Shape(int id) const; + std::unique_ptr Region() const; + void GetCellUnionBound(std::vector* cell_ids) const; + + const std::unique_ptr& Polygon() const { + return polygon_; + } + +private: + std::unique_ptr polygon_; +}; + + +// An S2Geography wrapping zero or more S2Geography objects. These objects +// can be used to represent a simple features GEOMETRYCOLLECTION. +class S2GeographyCollection: public S2Geography { +public: + S2GeographyCollection(): total_shapes_(0) {} + + S2GeographyCollection(std::vector> features): + features_(std::move(features)), + total_shapes_(0) { + for (const auto& feature: features_) { + num_shapes_.push_back(feature->num_shapes()); + total_shapes_ += feature->num_shapes(); + } + } + + int num_shapes() const; + std::unique_ptr Shape(int id) const; + std::unique_ptr Region() const; + + const std::vector>& Features() const { + return features_; + } + +private: + std::vector> features_; + std::vector num_shapes_; + int total_shapes_; +}; + + +// An S2Geography with a MutableS2ShapeIndex as the underlying data. +// These are used as inputs for operations that are implemented in S2 +// using the S2ShapeIndex (e.g., boolean operations). If an S2Geography +// instance will be used repeatedly, it will be faster to construct +// one ShapeIndexGeography and use it repeatedly. This class does not +// own any S2Geography objects that are added do it and thus is only +// valid for the scope of those objects. +class ShapeIndexGeography: public S2Geography { +public: + ShapeIndexGeography(MutableS2ShapeIndex::Options options = MutableS2ShapeIndex::Options()) + : shape_index_(options) {} + + explicit ShapeIndexGeography(const S2Geography& geog) { + Add(geog); + } + + // Add a S2Geography to the index, returning the last shape_id + // that was added to the index or -1 if no shapes were added + // to the index. + int Add(const S2Geography& geog) { + int id = -1; + for (int i = 0; i < geog.num_shapes(); i++) { + id = shape_index_.Add(geog.Shape(i)); + } + return id; + } + + int num_shapes() const; + std::unique_ptr Shape(int id) const; + std::unique_ptr Region() const; + + const MutableS2ShapeIndex& ShapeIndex() const { + return shape_index_; + } + +private: + MutableS2ShapeIndex shape_index_; +}; + + +} diff --git a/src/s2-geography/index.hpp b/src/s2-geography/index.hpp new file mode 100644 index 00000000..6715d8d0 --- /dev/null +++ b/src/s2-geography/index.hpp @@ -0,0 +1,96 @@ + +#pragma once + +#include + +#include "geography.hpp" + +namespace s2geography { + +// Unlike the ShapeIndexGeography, whose function is to index a single S2Geography +// or index multiple S2Geography objects as if they were a single S2Geography, +// the GeographyIndex exists to index a vector of S2Geography objects (like a +// GEOSSTRTree index), providing (hopefully) rapid access to possibly intersecting +// features. +class GeographyIndex { +public: + GeographyIndex(MutableS2ShapeIndex::Options options = MutableS2ShapeIndex::Options()) + : index_(options) {} + + void Add(const S2Geography& geog, int value) { + values_.reserve(values_.size() + geog.num_shapes()); + for (int i = 0; i < geog.num_shapes(); i++) { + int new_shape_id = index_.Add(geog.Shape(i)); + values_.resize(new_shape_id + 1); + values_[new_shape_id] = value; + } + } + + int value(int shape_id) const { + return values_[shape_id]; + } + + const MutableS2ShapeIndex& ShapeIndex() const { + return index_; + } + + MutableS2ShapeIndex& MutableShapeIndex() { + return index_; + } + + class Iterator { + public: + Iterator(const GeographyIndex* index): + index_(index), iterator_(&index_->ShapeIndex()) {} + + void Query(const std::vector& covering, std::unordered_set* indices) { + for (const S2CellId& query_cell: covering) { + Query(query_cell, indices); + } + } + + void Query(const S2CellId& cell_id, std::unordered_set* indices) { + S2ShapeIndex::CellRelation relation = iterator_.Locate(cell_id); + + if (relation == S2ShapeIndex::CellRelation::INDEXED) { + // We're in luck! these indexes have this cell in common + // add all the shapes it contains as possible intersectors + const S2ShapeIndexCell& index_cell = iterator_.cell(); + for (int k = 0; k < index_cell.num_clipped(); k++) { + int shape_id = index_cell.clipped(k).shape_id(); + indices->insert(index_->value(shape_id)); + } + } else if(relation == S2ShapeIndex::CellRelation::SUBDIVIDED) { + // Promising! the index has a child cell of iterator_.id() + // (at which iterator_ is now positioned). Keep iterating until the + // iterator is done OR we're no longer at a child cell of + // iterator_.id(). The ordering of the iterator isn't guaranteed anywhere + // in the documentation; however, this ordering would be consistent + // with that of a Normalized S2CellUnion. + while (!iterator_.done() && cell_id.contains(iterator_.id())) { + // add all the shapes the child cell contains as possible intersectors + const S2ShapeIndexCell& index_cell = iterator_.cell(); + for (int k = 0; k < index_cell.num_clipped(); k++) { + int shape_id = index_cell.clipped(k).shape_id(); + indices->insert(index_->value(shape_id)); + } + + // go to the next cell in the index + iterator_.Next(); + } + } + + // else: relation == S2ShapeIndex::CellRelation::DISJOINT (do nothing) + } + + private: + const GeographyIndex* index_; + MutableS2ShapeIndex::Iterator iterator_; + }; + +private: + MutableS2ShapeIndex index_; + std::vector values_; +}; + +} diff --git a/src/s2-geography/linear-referencing.cpp b/src/s2-geography/linear-referencing.cpp new file mode 100644 index 00000000..09f70bac --- /dev/null +++ b/src/s2-geography/linear-referencing.cpp @@ -0,0 +1,74 @@ + +#include "geography.hpp" +#include "linear-referencing.hpp" +#include "build.hpp" +#include "accessors.hpp" + +namespace s2geography { + +double s2_project_normalized(const PolylineGeography& geog1, + const S2Point& point) { + if (geog1.Polylines().size() != 1 || point.Norm2() == 0) { + return NAN; + } + + int next_vertex; + S2Point point_on_line = geog1.Polylines()[0]->Project(point, &next_vertex); + return geog1.Polylines()[0]->UnInterpolate(point_on_line, next_vertex); +} + +double s2_project_normalized(const S2Geography& geog1, const S2Geography& geog2) { + if (geog1.dimension() != 1 || geog2.dimension() != 0) { + return NAN; + } + + S2Point point; + for (int i = 0; i < geog2.num_shapes(); i++) { + auto shape = geog2.Shape(i); + for (int j = 0; j < shape->num_edges(); j++) { + if (point.Norm2() != 0) { + return NAN; + } else { + point = shape->edge(j).v0; + } + } + } + + auto geog1_poly_ptr = dynamic_cast(&geog1); + if (geog1_poly_ptr != nullptr) { + return s2_project_normalized(*geog1_poly_ptr, point); + } + + std::unique_ptr geog_poly = s2_rebuild(geog1, GlobalOptions()); + return s2_project_normalized(*geog_poly, geog2); +} + +S2Point s2_interpolate_normalized(const PolylineGeography& geog, double distance_norm) { + if (s2_is_empty(geog)) { + return S2Point(); + } else if (geog.Polylines().size() == 1) { + return geog.Polylines()[0]->Interpolate(distance_norm); + } else { + throw Exception("`geog` must contain 0 or 1 polyines"); + } +} + +S2Point s2_interpolate_normalized(const S2Geography& geog, double distance_norm) { + if (s2_is_empty(geog)) { + return S2Point(); + } + + if (geog.dimension() != 1 || geog.num_shapes() > 1) { + throw Exception("`geog` must be a single polyline"); + } + + auto geog_poly_ptr = dynamic_cast(&geog); + if (geog_poly_ptr != nullptr) { + return s2_interpolate_normalized(*geog_poly_ptr, distance_norm); + } + + std::unique_ptr geog_poly = s2_rebuild(geog, GlobalOptions()); + return s2_interpolate_normalized(*geog_poly, distance_norm); +} + +} diff --git a/src/s2-geography/linear-referencing.hpp b/src/s2-geography/linear-referencing.hpp new file mode 100644 index 00000000..7dba4a1e --- /dev/null +++ b/src/s2-geography/linear-referencing.hpp @@ -0,0 +1,11 @@ + +#pragma once + +#include "geography.hpp" + +namespace s2geography { + +double s2_project_normalized(const S2Geography& geog1, const S2Geography& geog2); +S2Point s2_interpolate_normalized(const S2Geography& geog, double distance_norm); + +} diff --git a/src/s2-geography/predicates.cpp b/src/s2-geography/predicates.cpp new file mode 100644 index 00000000..f2003f9c --- /dev/null +++ b/src/s2-geography/predicates.cpp @@ -0,0 +1,95 @@ + +#include "s2/s2boolean_operation.h" +#include "s2/s2edge_tessellator.h" +#include "s2/s2lax_loop_shape.h" + +#include "accessors.hpp" +#include "predicates.hpp" + +namespace s2geography { + +bool s2_intersects(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + const S2BooleanOperation::Options& options) { + return S2BooleanOperation::Intersects( + geog1.ShapeIndex(), + geog2.ShapeIndex(), + options); +} + +bool s2_equals(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + const S2BooleanOperation::Options& options) { + return S2BooleanOperation::Equals( + geog1.ShapeIndex(), + geog2.ShapeIndex(), + options); +} + +bool s2_contains(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + const S2BooleanOperation::Options& options) { + if (s2_is_empty(geog2)) { + return false; + } else { + return S2BooleanOperation::Contains( + geog1.ShapeIndex(), + geog2.ShapeIndex(), + options); + } +} + +// Note that 'touches' can be implemeted using: +// +// S2BooleanOperation::Options closedOptions = options; +// closedOptions.set_polygon_model(S2BooleanOperation::PolygonModel::CLOSED); +// closedOptions.set_polyline_model(S2BooleanOperation::PolylineModel::CLOSED); +// S2BooleanOperation::Options openOptions = options; +// openOptions.set_polygon_model(S2BooleanOperation::PolygonModel::OPEN); +// openOptions.set_polyline_model(S2BooleanOperation::PolylineModel::OPEN); +// s2_intersects(geog1, geog2, closed_options) && +// !s2_intersects(geog1, geog2, open_options); +// +// ...it isn't implemented here because the options creation should be done +// outside of any loop. + +bool s2_intersects_box(const ShapeIndexGeography& geog1, + const S2LatLngRect& rect, + const S2BooleanOperation::Options& options, + double tolerance) { + // 99% of this is making a S2Loop out of a S2LatLngRect + // This should probably be implemented elsewhere + S2::PlateCarreeProjection projection(180); + S2EdgeTessellator tessellator(&projection, S1Angle::Degrees(tolerance)); + std::vector vertices; + + tessellator.AppendUnprojected( + R2Point(rect.lng_lo().degrees(), rect.lat_lo().degrees()), + R2Point(rect.lng_hi().degrees(), rect.lat_lo().degrees()), + &vertices); + tessellator.AppendUnprojected( + R2Point(rect.lng_hi().degrees(), rect.lat_lo().degrees()), + R2Point(rect.lng_hi().degrees(), rect.lat_hi().degrees()), + &vertices); + tessellator.AppendUnprojected( + R2Point(rect.lng_hi().degrees(), rect.lat_hi().degrees()), + R2Point(rect.lng_lo().degrees(), rect.lat_hi().degrees()), + &vertices); + tessellator.AppendUnprojected( + R2Point(rect.lng_lo().degrees(), rect.lat_hi().degrees()), + R2Point(rect.lng_lo().degrees(), rect.lat_lo().degrees()), + &vertices); + + vertices.pop_back(); + + auto loop = absl::make_unique(std::move(vertices)); + MutableS2ShapeIndex index; + index.Add(std::move(loop)); + + return S2BooleanOperation::Intersects( + geog1.ShapeIndex(), + index, + options); +} + +} diff --git a/src/s2-geography/predicates.hpp b/src/s2-geography/predicates.hpp new file mode 100644 index 00000000..c267d733 --- /dev/null +++ b/src/s2-geography/predicates.hpp @@ -0,0 +1,31 @@ + +#pragma once + +#include "s2/s2boolean_operation.h" + +#include "geography.hpp" + +namespace s2geography { + +bool s2_intersects(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + const S2BooleanOperation::Options& options); + +bool s2_equals(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + const S2BooleanOperation::Options& options); + +bool s2_contains(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + const S2BooleanOperation::Options& options); + +bool s2_touches(const ShapeIndexGeography& geog1, + const ShapeIndexGeography& geog2, + const S2BooleanOperation::Options& options); + +bool s2_intersects_box(const ShapeIndexGeography& geog1, + const S2LatLngRect& rect, + const S2BooleanOperation::Options& options, + double tolerance); + +} diff --git a/src/s2-geography/s2-geography.hpp b/src/s2-geography/s2-geography.hpp new file mode 100644 index 00000000..39793a20 --- /dev/null +++ b/src/s2-geography/s2-geography.hpp @@ -0,0 +1,13 @@ + +#pragma once + +#include "geography.hpp" +#include "index.hpp" +#include "accessors.hpp" +#include "accessors-geog.hpp" +#include "linear-referencing.hpp" +#include "distance.hpp" +#include "build.hpp" +#include "coverings.hpp" +#include "predicates.hpp" +#include "constructor.hpp" diff --git a/src/s2-lnglat.cpp b/src/s2-lnglat.cpp index 764e5ccd..b5153a10 100644 --- a/src/s2-lnglat.cpp +++ b/src/s2-lnglat.cpp @@ -1,9 +1,6 @@ #include "s2/s2latlng.h" #include "s2/s2point.h" -#include "wk/rcpp-io.hpp" -#include "wk/wkb-reader.hpp" -#include "wk/wkb-writer.hpp" #include using namespace Rcpp; diff --git a/src/s2-matrix.cpp b/src/s2-matrix.cpp index 10a79735..84c13dad 100644 --- a/src/s2-matrix.cpp +++ b/src/s2-matrix.cpp @@ -14,106 +14,41 @@ #include using namespace Rcpp; -std::unordered_map buildSourcedIndex(List geog, MutableS2ShapeIndex* index) { - std::unordered_map indexSource; - std::vector shapeIds; - - for (R_xlen_t j = 0; j < geog.size(); j++) { - checkUserInterrupt(); - SEXP item2 = geog[j]; - - // build index and store index IDs so that shapeIds can be - // mapped back to the geog index - if (item2 == R_NilValue) { - Rcpp::stop("Missing `y` not allowed in binary indexed operators()"); - } else { - Rcpp::XPtr feature2(item2); - shapeIds = feature2->BuildShapeIndex(index); - for (size_t k = 0; k < shapeIds.size(); k ++) { - indexSource[shapeIds[k]] = j; - } - } - } - - return indexSource; -} - -std::unordered_set findPossibleIntersections(const S2Region& region, - const MutableS2ShapeIndex* index, - std::unordered_map& source, - int maxRegionCells) { - - std::unordered_set mightIntersectIndices; - MutableS2ShapeIndex::Iterator indexIterator(index); - - // generate a small covering of the region - S2RegionCoverer coverer; - coverer.mutable_options()->set_max_cells(maxRegionCells); - S2CellUnion covering = coverer.GetCovering(region); - - // iterate over cells in the featureIndex - for (S2CellId featureCellId: covering) { - S2ShapeIndex::CellRelation relation = indexIterator.Locate(featureCellId); - - if (relation == S2ShapeIndex::CellRelation::INDEXED) { - // we're in luck! these indexes have this cell in common - // add all the features it contains as possible intersectors for featureIndex - const S2ShapeIndexCell& cell = indexIterator.cell(); - for (int k = 0; k < cell.num_clipped(); k++) { - int shapeId = cell.clipped(k).shape_id(); - mightIntersectIndices.insert(source[shapeId]); - } - - } else if(relation == S2ShapeIndex::CellRelation::SUBDIVIDED) { - // promising! the geog2 index has a child cell of it.id() - // (at which indexIterator is now positioned) - // keep iterating until the iterator is done OR we're no longer at a child cell of - // it.id(). The ordering of the iterator isn't guaranteed anywhere in the documentation; - // however, this ordering would be consistent with that of a Normalized - // S2CellUnion. - while (!indexIterator.done() && featureCellId.contains(indexIterator.id())) { - // potentially many cells in the indexIterator, so let the user cancel if this is - // running too long - checkUserInterrupt(); - - // add all the features the child cell contains as possible intersectors for featureIndex - const S2ShapeIndexCell& cell = indexIterator.cell(); - for (int k = 0; k < cell.num_clipped(); k++) { - int shapeId = cell.clipped(k).shape_id(); - mightIntersectIndices.insert(source[shapeId]); - } - - // go to the next cell in the index - indexIterator.Next(); - } - } - - // else: relation == S2ShapeIndex::CellRelation::DISJOINT (do nothing) - } - - return mightIntersectIndices; -} template class IndexedBinaryGeographyOperator: public UnaryGeographyOperator { public: - std::unique_ptr geog2Index; - std::unordered_map geog2IndexSource; - - IndexedBinaryGeographyOperator() { - this->geog2Index = absl::make_unique(); - } + std::unique_ptr geog2_index; + std::unique_ptr iterator; - // maxEdgesPerCell should be between 10 and 50, with lower numbers + // max_edges_per_cell should be between 10 and 50, with lower numbers // leading to more memory usage (but potentially faster query times). Benchmarking // with binary prediates seems to indicate that values on the high end // of the spectrum do a reasonable job of efficient preselection, and that // decreasing this value does little to increase performance. - virtual void buildIndex(List geog2, int maxEdgesPerCell = 50) { - MutableS2ShapeIndex::Options indexOptions; - indexOptions.set_max_edges_per_cell(maxEdgesPerCell); - this->geog2Index = absl::make_unique(indexOptions); - this->geog2IndexSource = buildSourcedIndex(geog2, this->geog2Index.get()); + + IndexedBinaryGeographyOperator(int maxEdgesPerCell = 50) { + MutableS2ShapeIndex::Options index_options; + index_options.set_max_edges_per_cell(maxEdgesPerCell); + geog2_index = absl::make_unique(index_options); + } + + virtual void buildIndex(List geog2) { + for (R_xlen_t j = 0; j < geog2.size(); j++) { + checkUserInterrupt(); + SEXP item2 = geog2[j]; + + // build index and store index IDs so that shapeIds can be + // mapped back to the geog index + if (item2 == R_NilValue) { + Rcpp::stop("Missing `y` not allowed in binary indexed operators()"); + } else { + Rcpp::XPtr feature2(item2); + geog2_index->Add(feature2->Geog(), j); + } + } + + iterator = absl::make_unique(geog2_index.get()); } }; @@ -125,14 +60,14 @@ IntegerVector cpp_s2_closest_feature(List geog1, List geog2) { class Op: public IndexedBinaryGeographyOperator { public: int processFeature(Rcpp::XPtr feature, R_xlen_t i) { - S2ClosestEdgeQuery query(this->geog2Index.get()); - S2ClosestEdgeQuery::ShapeIndexTarget target(feature->ShapeIndex()); + S2ClosestEdgeQuery query(&geog2_index->ShapeIndex()); + S2ClosestEdgeQuery::ShapeIndexTarget target(&feature->Index().ShapeIndex()); const auto& result = query.FindClosestEdge(&target); if (result.is_empty()) { return NA_INTEGER; } else { // convert to R index (+1) - return this->geog2IndexSource[result.shape_id()] + 1; + return geog2_index->value(result.shape_id()) + 1; } } }; @@ -148,14 +83,14 @@ IntegerVector cpp_s2_farthest_feature(List geog1, List geog2) { class Op: public IndexedBinaryGeographyOperator { public: int processFeature(Rcpp::XPtr feature, R_xlen_t i) { - S2FurthestEdgeQuery query(this->geog2Index.get()); - S2FurthestEdgeQuery::ShapeIndexTarget target(feature->ShapeIndex()); + S2FurthestEdgeQuery query(&geog2_index->ShapeIndex()); + S2FurthestEdgeQuery::ShapeIndexTarget target(&feature->Index().ShapeIndex()); const auto& result = query.FindFurthestEdge(&target); if (result.is_empty()) { return NA_INTEGER; } else { // convert to R index (+1) - return this->geog2IndexSource[result.shape_id()] + 1; + return geog2_index->value(result.shape_id()) + 1; } } }; @@ -172,17 +107,17 @@ List cpp_s2_closest_edges(List geog1, List geog2, int n, double min_distance, class Op: public IndexedBinaryGeographyOperator { public: IntegerVector processFeature(Rcpp::XPtr feature, R_xlen_t i) { - S2ClosestEdgeQuery query(this->geog2Index.get()); + S2ClosestEdgeQuery query(&geog2_index->ShapeIndex()); query.mutable_options()->set_max_results(n); query.mutable_options()->set_max_distance(S1ChordAngle::Radians(max_distance)); - S2ClosestEdgeQuery::ShapeIndexTarget target(feature->ShapeIndex()); + S2ClosestEdgeQuery::ShapeIndexTarget target(&feature->Index().ShapeIndex()); const auto& result = query.FindClosestEdges(&target); // this code searches edges, which may come from the same feature std::unordered_set features; for (S2ClosestEdgeQuery::Result res : result) { if (res.distance().radians() > this->min_distance) { - features.insert(this->geog2IndexSource[res.shape_id()] + 1); + features.insert(geog2_index->value(res.shape_id()) + 1); } } @@ -210,56 +145,57 @@ class IndexedMatrixPredicateOperator: public IndexedBinaryGeographyOperator(maxEdgesPerCell), maxFeatureCells(maxFeatureCells) { GeographyOperationOptions options(s2options); this->options = options.booleanOperationOptions(); + this->coverer.mutable_options()->set_max_cells(maxFeatureCells); } - // See IndexedBinaryGeographyOperator::buildIndex() for why 50 is the default value - // for maxEdgesPerCell - void buildIndex(List geog2, int maxEdgesPerCell = 50) { + void buildIndex(List geog2) { this->geog2 = geog2; - IndexedBinaryGeographyOperator::buildIndex(geog2, maxEdgesPerCell); + IndexedBinaryGeographyOperator::buildIndex(geog2); } IntegerVector processFeature(Rcpp::XPtr feature, R_xlen_t i) { - S2ShapeIndex* index1 = feature->ShapeIndex(); - S2ShapeIndexRegion region = MakeS2ShapeIndexRegion(index1); - - // build a list of candidate feature indices - std::unordered_set mightIntersectIndices = findPossibleIntersections( - region, - this->geog2Index.get(), - this->geog2IndexSource, - this->maxFeatureCells - ); + coverer.GetCovering(*feature->Geog().Region(), &cell_ids); + indices_unsorted.clear(); + iterator->Query(cell_ids, &indices_unsorted); // loop through features from geog2 that might intersect feature // and build a list of indices that actually intersect (based on // this->actuallyIntersects(), which might perform alternative // comparisons) - std::vector actuallyIntersectIndices; - for (R_xlen_t j: mightIntersectIndices) { + indices.clear(); + for (int j: indices_unsorted) { SEXP item = this->geog2[j]; XPtr feature2(item); - if (this->actuallyIntersects(index1, feature2->ShapeIndex(), i, j)) { + + if (this->actuallyIntersects(feature->Index(), feature2->Index(), i, j)) { // convert to R index here + 1 - actuallyIntersectIndices.push_back(j + 1); + indices.push_back(j + 1); } } // return sorted integer vector - std::sort(actuallyIntersectIndices.begin(), actuallyIntersectIndices.end()); - return Rcpp::IntegerVector(actuallyIntersectIndices.begin(), actuallyIntersectIndices.end()); + std::sort(indices.begin(), indices.end()); + return Rcpp::IntegerVector(indices.begin(), indices.end()); }; - virtual bool actuallyIntersects(S2ShapeIndex* index1, S2ShapeIndex* index2, R_xlen_t i, R_xlen_t j) = 0; + virtual bool actuallyIntersects(const s2geography::ShapeIndexGeography& index1, + const s2geography::ShapeIndexGeography& index2, + R_xlen_t i, R_xlen_t j) = 0; protected: List geog2; S2BooleanOperation::Options options; int maxFeatureCells; + S2RegionCoverer coverer; + std::vector cell_ids; + std::unordered_set indices_unsorted; + std::vector indices; }; // [[Rcpp::export]] @@ -267,16 +203,18 @@ List cpp_s2_may_intersect_matrix(List geog1, List geog2, int maxEdgesPerCell, int maxFeatureCells, List s2options) { class Op: public IndexedMatrixPredicateOperator { public: - Op(List s2options, int maxFeatureCells): - IndexedMatrixPredicateOperator(s2options, maxFeatureCells) {} + Op(List s2options, int maxFeatureCells, int maxEdgesPerCell): + IndexedMatrixPredicateOperator(s2options, maxFeatureCells, maxEdgesPerCell) {} - bool actuallyIntersects(S2ShapeIndex* index1, S2ShapeIndex* index2, R_xlen_t i, R_xlen_t j) { + bool actuallyIntersects(const s2geography::ShapeIndexGeography& index1, + const s2geography::ShapeIndexGeography& index2, + R_xlen_t i, R_xlen_t j) { return true; }; }; - Op op(s2options, maxFeatureCells); - op.buildIndex(geog2, maxEdgesPerCell); + Op op(s2options, maxFeatureCells, maxEdgesPerCell); + op.buildIndex(geog2); return op.processVector(geog1); } @@ -285,8 +223,10 @@ List cpp_s2_contains_matrix(List geog1, List geog2, List s2options) { class Op: public IndexedMatrixPredicateOperator { public: Op(List s2options): IndexedMatrixPredicateOperator(s2options) {} - bool actuallyIntersects(S2ShapeIndex* index1, S2ShapeIndex* index2, R_xlen_t i, R_xlen_t j) { - return S2BooleanOperation::Contains(*index1, *index2, this->options); + bool actuallyIntersects(const s2geography::ShapeIndexGeography& index1, + const s2geography::ShapeIndexGeography& index2, + R_xlen_t i, R_xlen_t j) { + return s2geography::s2_contains(index1, index2, this->options); }; }; @@ -300,9 +240,11 @@ List cpp_s2_within_matrix(List geog1, List geog2, List s2options) { class Op: public IndexedMatrixPredicateOperator { public: Op(List s2options): IndexedMatrixPredicateOperator(s2options) {} - bool actuallyIntersects(S2ShapeIndex* index1, S2ShapeIndex* index2, R_xlen_t i, R_xlen_t j) { + bool actuallyIntersects(const s2geography::ShapeIndexGeography& index1, + const s2geography::ShapeIndexGeography& index2, + R_xlen_t i, R_xlen_t j) { // note reversed index2, index1 - return S2BooleanOperation::Contains(*index2, *index1, this->options); + return s2geography::s2_contains(index2, index1, this->options); }; }; @@ -316,8 +258,10 @@ List cpp_s2_intersects_matrix(List geog1, List geog2, List s2options) { class Op: public IndexedMatrixPredicateOperator { public: Op(List s2options): IndexedMatrixPredicateOperator(s2options) {} - bool actuallyIntersects(S2ShapeIndex* index1, S2ShapeIndex* index2, R_xlen_t i, R_xlen_t j) { - return S2BooleanOperation::Intersects(*index1, *index2, this->options); + bool actuallyIntersects(const s2geography::ShapeIndexGeography& index1, + const s2geography::ShapeIndexGeography& index2, + R_xlen_t i, R_xlen_t j) { + return s2geography::s2_intersects(index1, index2, this->options); }; }; @@ -331,8 +275,10 @@ List cpp_s2_equals_matrix(List geog1, List geog2, List s2options) { class Op: public IndexedMatrixPredicateOperator { public: Op(List s2options): IndexedMatrixPredicateOperator(s2options) {} - bool actuallyIntersects(S2ShapeIndex* index1, S2ShapeIndex* index2, R_xlen_t i, R_xlen_t j) { - return S2BooleanOperation::Equals(*index1, *index2, this->options); + bool actuallyIntersects(const s2geography::ShapeIndexGeography& index1, + const s2geography::ShapeIndexGeography& index2, + R_xlen_t i, R_xlen_t j) { + return s2geography::s2_equals(index1, index2, this->options); }; }; @@ -355,10 +301,11 @@ List cpp_s2_touches_matrix(List geog1, List geog2, List s2options) { this->openOptions.set_polyline_model(S2BooleanOperation::PolylineModel::OPEN); } - bool actuallyIntersects(S2ShapeIndex* index1, S2ShapeIndex* index2, R_xlen_t i, R_xlen_t j) { - // efficiently re-uses the index on geog2 and takes advantage of short-circuiting && - return S2BooleanOperation::Intersects(*index1, *index2, this->closedOptions) && - !S2BooleanOperation::Intersects(*index1, *index2, this->openOptions); + bool actuallyIntersects(const s2geography::ShapeIndexGeography& index1, + const s2geography::ShapeIndexGeography& index2, + R_xlen_t i, R_xlen_t j) { + return s2geography::s2_intersects(index1, index2, this->closedOptions) && + !s2geography::s2_intersects(index1, index2, this->openOptions); }; private: @@ -441,8 +388,8 @@ List cpp_s2_dwithin_matrix(List geog1, List geog2, double distance) { Op(double distance): distance(distance) {} bool processFeature(XPtr feature1, XPtr feature2, R_xlen_t i, R_xlen_t j) { - S2ClosestEdgeQuery query(feature2->ShapeIndex()); - S2ClosestEdgeQuery::ShapeIndexTarget target(feature1->ShapeIndex()); + S2ClosestEdgeQuery query(&feature2->Index().ShapeIndex()); + S2ClosestEdgeQuery::ShapeIndexTarget target(&feature1->Index().ShapeIndex()); return query.IsDistanceLessOrEqual(&target, S1ChordAngle::Radians(this->distance)); }; }; @@ -500,8 +447,8 @@ NumericMatrix cpp_s2_distance_matrix(List geog1, List geog2) { double processFeature(XPtr feature1, XPtr feature2, R_xlen_t i, R_xlen_t j) { - S2ClosestEdgeQuery query(feature1->ShapeIndex()); - S2ClosestEdgeQuery::ShapeIndexTarget target(feature2->ShapeIndex()); + S2ClosestEdgeQuery query(&feature1->Index().ShapeIndex()); + S2ClosestEdgeQuery::ShapeIndexTarget target(&feature2->Index().ShapeIndex()); const auto& result = query.FindClosestEdge(&target); S1ChordAngle angle = result.distance(); @@ -525,8 +472,8 @@ NumericMatrix cpp_s2_max_distance_matrix(List geog1, List geog2) { double processFeature(XPtr feature1, XPtr feature2, R_xlen_t i, R_xlen_t j) { - S2FurthestEdgeQuery query(feature1->ShapeIndex()); - S2FurthestEdgeQuery::ShapeIndexTarget target(feature2->ShapeIndex()); + S2FurthestEdgeQuery query(&feature1->Index().ShapeIndex()); + S2FurthestEdgeQuery::ShapeIndexTarget target(&feature2->Index().ShapeIndex()); const auto& result = query.FindFurthestEdge(&target); S1ChordAngle angle = result.distance(); @@ -557,17 +504,7 @@ List cpp_s2_contains_matrix_brute_force(List geog1, List geog2, List s2options) Op(List s2options): BruteForceMatrixPredicateOperator(s2options) {} bool processFeature(XPtr feature1, XPtr feature2, R_xlen_t i, R_xlen_t j) { - // by default Contains() will return true for Contains(x, EMPTY), which is - // not true in BigQuery or GEOS - if (feature2->IsEmpty()) { - return false; - } else { - return S2BooleanOperation::Contains( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - this->options - ); - } + return s2geography::s2_contains(feature1->Index(), feature2->Index(), options); }; }; @@ -583,18 +520,7 @@ List cpp_s2_within_matrix_brute_force(List geog1, List geog2, List s2options) { bool processFeature(XPtr feature1, XPtr feature2, R_xlen_t i, R_xlen_t j) { // note reversed index2, index1 - - // by default Contains() will return true for Contains(x, EMPTY), which is - // not true in BigQuery or GEOS - if (feature1->IsEmpty()) { - return false; - } else { - return S2BooleanOperation::Contains( - *feature2->ShapeIndex(), - *feature1->ShapeIndex(), - this->options - ); - } + return s2geography::s2_contains(feature2->Index(), feature1->Index(), options); }; }; @@ -609,11 +535,7 @@ List cpp_s2_intersects_matrix_brute_force(List geog1, List geog2, List s2options Op(List s2options): BruteForceMatrixPredicateOperator(s2options) {} bool processFeature(XPtr feature1, XPtr feature2, R_xlen_t i, R_xlen_t j) { - return S2BooleanOperation::Intersects( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - this->options - ); + return s2geography::s2_intersects(feature1->Index(), feature2->Index(), options); } }; @@ -628,11 +550,7 @@ List cpp_s2_disjoint_matrix_brute_force(List geog1, List geog2, List s2options) Op(List s2options): BruteForceMatrixPredicateOperator(s2options) {} bool processFeature(XPtr feature1, XPtr feature2, R_xlen_t i, R_xlen_t j) { - return !S2BooleanOperation::Intersects( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - this->options - ); + return !s2geography::s2_intersects(feature1->Index(), feature2->Index(), options); } }; @@ -647,11 +565,7 @@ List cpp_s2_equals_matrix_brute_force(List geog1, List geog2, List s2options) { Op(List s2options): BruteForceMatrixPredicateOperator(s2options) {} bool processFeature(XPtr feature1, XPtr feature2, R_xlen_t i, R_xlen_t j) { - return S2BooleanOperation::Equals( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - this->options - ); + return s2geography::s2_equals(feature1->Index(), feature2->Index(), options); } }; diff --git a/src/s2-options.h b/src/s2-options.h index 778628d8..86ec1c18 100644 --- a/src/s2-options.h +++ b/src/s2-options.h @@ -10,6 +10,8 @@ #include "s2/s2builderutil_s2polyline_vector_layer.h" #include "s2/s2builderutil_s2point_vector_layer.h" +#include "s2-geography/s2-geography.hpp" + // This class wraps several concepts in the S2BooleanOperation, // and S2Layer, parameterized such that these can be specified from R class GeographyOperationOptions { @@ -194,6 +196,32 @@ class GeographyOperationOptions { return options; } + // options for new GlobalOptions API + s2geography::GlobalOptions geographyOptions() { + s2geography::GlobalOptions options; + options.boolean_operation = booleanOperationOptions(); + options.builder = builderOptions(); + + LayerOptions layer_options = layerOptions(); + options.point_layer = layer_options.pointLayerOptions; + options.polyline_layer = layer_options.polylineLayerOptions; + options.polygon_layer = layer_options.polygonLayerOptions; + + if (!(layer_options.dimensions & Dimension::POINT)) { + options.point_layer_action = s2geography::GlobalOptions::OUTPUT_ACTION_IGNORE; + } + + if (!(layer_options.dimensions & Dimension::POLYLINE)) { + options.polyline_layer_action = s2geography::GlobalOptions::OUTPUT_ACTION_IGNORE; + } + + if (!(layer_options.dimensions & Dimension::POLYGON)) { + options.polygon_layer_action = s2geography::GlobalOptions::OUTPUT_ACTION_IGNORE; + } + + return options; + } + // build options for S2Builder S2Builder::Options builderOptions() { S2Builder::Options options; diff --git a/src/s2-predicates.cpp b/src/s2-predicates.cpp index 67aea0bf..3134b5bb 100644 --- a/src/s2-predicates.cpp +++ b/src/s2-predicates.cpp @@ -28,11 +28,7 @@ LogicalVector cpp_s2_intersects(List geog1, List geog2, List s2options) { public: Op(List s2options): BinaryPredicateOperator(s2options) {} int processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - return S2BooleanOperation::Intersects( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - options - ); + return s2geography::s2_intersects(feature1->Index(), feature2->Index(), options); }; }; @@ -47,11 +43,7 @@ LogicalVector cpp_s2_equals(List geog1, List geog2, List s2options) { public: Op(List s2options): BinaryPredicateOperator(s2options) {} int processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - return S2BooleanOperation::Equals( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - this->options - ); + return s2geography::s2_equals(feature1->Index(), feature2->Index(), options); } }; @@ -65,17 +57,7 @@ LogicalVector cpp_s2_contains(List geog1, List geog2, List s2options) { public: Op(List s2options): BinaryPredicateOperator(s2options) {} int processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - // by default Contains() will return true for Contains(x, EMPTY), which is - // not true in BigQuery or GEOS - if (feature2->IsEmpty()) { - return false; - } else { - return S2BooleanOperation::Contains( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - this->options - ); - } + return s2geography::s2_contains(feature1->Index(), feature2->Index(), options); } }; @@ -98,16 +80,8 @@ LogicalVector cpp_s2_touches(List geog1, List geog2, List s2options) { } int processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - return S2BooleanOperation::Intersects( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - this->closedOptions - ) && - !S2BooleanOperation::Intersects( - *feature1->ShapeIndex(), - *feature2->ShapeIndex(), - this->openOptions - ); + return s2geography::s2_intersects(feature1->Index(), feature2->Index(), this->closedOptions) && + !s2geography::s2_intersects(feature1->Index(), feature2->Index(), this->openOptions); } private: @@ -131,8 +105,8 @@ LogicalVector cpp_s2_dwithin(List geog1, List geog2, NumericVector distance) { Op(NumericVector distance): distance(distance) {} int processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - S2ClosestEdgeQuery query(feature1->ShapeIndex()); - S2ClosestEdgeQuery::ShapeIndexTarget target(feature2->ShapeIndex()); + S2ClosestEdgeQuery query(&feature1->Index().ShapeIndex()); + S2ClosestEdgeQuery::ShapeIndexTarget target(&feature2->Index().ShapeIndex()); return query.IsDistanceLessOrEqual(&target, S1ChordAngle::Radians(this->distance[i])); } }; @@ -190,39 +164,9 @@ LogicalVector cpp_s2_intersects_box(List geog, return false; } - // create polygon vertices - std::vector points(2 + 2 * detail); - S2LatLng vertex; - - // south edge - for (int i = 0; i <= detail; i++) { - vertex = S2LatLng::FromDegrees(xmin + deltaDegrees * i, ymin).Normalized(); - points[i] = vertex.ToPoint(); - } - - // north edge - for (int i = 0; i <= detail; i++) { - vertex = S2LatLng::FromDegrees(xmax - deltaDegrees * i, ymax).Normalized(); - points[detail + 1 + i] = vertex.ToPoint(); - } + S2LatLngRect rect(S2LatLng::FromDegrees(ymin, xmin), S2LatLng::FromDegrees(ymax, xmax)); - // create polygon - std::unique_ptr loop(new S2Loop()); - loop->set_s2debug_override(S2Debug::DISABLE); - loop->Init(points); - loop->Normalize(); - - std::vector> loops(1); - loops[0] = std::move(loop); - S2Polygon polygon; - polygon.InitOriented(std::move(loops)); - - // test intersection - return S2BooleanOperation::Intersects( - polygon.index(), - *feature->ShapeIndex(), - this->options - ); + return s2geography::s2_intersects_box(feature->Index(), rect, options, deltaDegrees); } }; diff --git a/src/s2-transformers.cpp b/src/s2-transformers.cpp index f09fbe9e..408d0b87 100644 --- a/src/s2-transformers.cpp +++ b/src/s2-transformers.cpp @@ -1,196 +1,34 @@ -#include "s2/s2boolean_operation.h" -#include "s2/s2closest_edge_query.h" -#include "s2/s2polygon.h" -#include "s2/s2polyline.h" -#include "s2/s2point.h" -#include "s2/s2error.h" -#include "s2/s2boolean_operation.h" -#include "s2/s2builder.h" -#include "s2/s2builderutil_s2polygon_layer.h" -#include "s2/s2builderutil_s2polyline_vector_layer.h" -#include "s2/s2builderutil_s2point_vector_layer.h" -#include "s2/s2builderutil_closed_set_normalizer.h" -#include "s2/s2builderutil_snap_functions.h" #include "s2/s2shape_index_buffered_region.h" #include "s2/s2region_coverer.h" -#include "s2/s2convex_hull_query.h" #include "s2-options.h" #include "geography-operator.h" -#include "point-geography.h" -#include "polyline-geography.h" -#include "polygon-geography.h" -#include "geography-collection.h" #include using namespace Rcpp; -std::unique_ptr geographyFromLayers(std::vector points, - std::vector> polylines, - std::unique_ptr polygon, - int dimensions) { - // count non-empty dimensions - bool has_polygon = (dimensions & GeographyOperationOptions::Dimension::POLYGON) && - !polygon->is_empty(); - bool has_polyline = (dimensions & GeographyOperationOptions::Dimension::POLYLINE) && - (polylines.size() > 0); - bool has_points = (dimensions & GeographyOperationOptions::Dimension::POINT) && - (points.size() > 0); - int nonEmptyDimensions = has_polygon + has_polyline + has_points; - - // return empty output - if (nonEmptyDimensions == 0) { - return absl::make_unique(); - } - - // return mixed dimension output - if (nonEmptyDimensions > 1) { - std::vector> features; - - if (has_points) { - features.push_back(absl::make_unique(std::move(points))); - } - - if (has_polyline) { - features.push_back(absl::make_unique(std::move(polylines))); - } - - if (has_polygon) { - features.push_back(absl::make_unique(std::move(polygon))); - } - - return absl::make_unique(std::move(features)); - } - - // return single dimension output - if (has_polygon) { - return absl::make_unique(std::move(polygon)); - } else if (has_polyline) { - return absl::make_unique(std::move(polylines)); - } else { - return absl::make_unique(std::move(points)); - } -} - -std::unique_ptr doBooleanOperation(S2ShapeIndex* index1, S2ShapeIndex* index2, - S2BooleanOperation::OpType opType, - S2BooleanOperation::Options options, - GeographyOperationOptions::LayerOptions layerOptions) { - - // create the data structures that will contain the output - std::vector points; - std::vector> polylines; - std::unique_ptr polygon = absl::make_unique(); - - s2builderutil::LayerVector layers(3); - layers[0] = absl::make_unique(&points, layerOptions.pointLayerOptions); - layers[1] = absl::make_unique(&polylines, layerOptions.polylineLayerOptions); - layers[2] = absl::make_unique(polygon.get(), layerOptions.polygonLayerOptions); - - // do the boolean operation - S2BooleanOperation booleanOp( - opType, - // normalizing the closed set here is required for line intersections - // to work as expected - s2builderutil::NormalizeClosedSet(std::move(layers)), - options - ); - - // build and check for errors - S2Error error; - if (!booleanOp.Build(*index1, *index2, &error)) { - stop(error.text()); - } - - // construct output - return geographyFromLayers( - std::move(points), - std::move(polylines), - std::move(polygon), - layerOptions.dimensions - ); -} - -std::unique_ptr rebuildGeography(S2ShapeIndex* index, - S2Builder::Options options, - GeographyOperationOptions::LayerOptions layerOptions) { - // create the builder - S2Builder builder(options); - - // create the data structures that will contain the output - std::vector points; - std::vector> polylines; - std::unique_ptr polygon = absl::make_unique(); - - // add shapes to the layer with the appropriate dimension - builder.StartLayer( - absl::make_unique(&points, layerOptions.pointLayerOptions) - ); - for (S2Shape* shape : *index) { - if (shape->dimension() == 0) { - builder.AddShape(*shape); - } - } - - builder.StartLayer( - absl::make_unique(&polylines, layerOptions.polylineLayerOptions) - ); - for (S2Shape* shape : *index) { - if (shape->dimension() == 1) { - builder.AddShape(*shape); - } - } - - builder.StartLayer( - absl::make_unique(polygon.get(), layerOptions.polygonLayerOptions) - ); - for (S2Shape* shape : *index) { - if (shape->dimension() == 2) { - builder.AddShape(*shape); - } - } - - // build the output - S2Error error; - if (!builder.Build(&error)) { - throw GeographyOperatorException(error.text()); - } - - // construct output - return geographyFromLayers( - std::move(points), - std::move(polylines), - std::move(polygon), - layerOptions.dimensions - ); -} class BooleanOperationOp: public BinaryGeographyOperator { public: BooleanOperationOp(S2BooleanOperation::OpType opType, List s2options): opType(opType) { GeographyOperationOptions options(s2options); - this->options = options.booleanOperationOptions(); - this->layerOptions = options.layerOptions(); + this->geography_options = options.geographyOptions(); } SEXP processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - std::unique_ptr geography = doBooleanOperation( - feature1->ShapeIndex(), - feature2->ShapeIndex(), + std::unique_ptr geog_out = s2geography::s2_boolean_operation( + feature1->Index(), feature2->Index(), this->opType, - this->options, - this->layerOptions - ); + this->geography_options); - return Rcpp::XPtr(geography.release()); + return Geography::MakeXPtr(std::move(geog_out)); } private: S2BooleanOperation::OpType opType; - S2BooleanOperation::Options options; - GeographyOperationOptions::LayerOptions layerOptions; + s2geography::GlobalOptions geography_options; }; // [[Rcpp::export]] @@ -220,8 +58,8 @@ List cpp_s2_sym_difference(List geog1, List geog2, List s2options) { // [[Rcpp::export]] List cpp_s2_coverage_union_agg(List geog, List s2options, bool naRm) { GeographyOperationOptions options(s2options); + s2geography::S2CoverageUnionAggregator agg(options.geographyOptions()); - MutableS2ShapeIndex index; SEXP item; for (R_xlen_t i = 0; i < geog.size(); i++) { item = geog[i]; @@ -231,36 +69,18 @@ List cpp_s2_coverage_union_agg(List geog, List s2options, bool naRm) { if (item != R_NilValue) { Rcpp::XPtr feature(item); - feature->BuildShapeIndex(&index); + agg.Add(feature->Geog()); } } - MutableS2ShapeIndex emptyIndex; - std::unique_ptr geography = doBooleanOperation( - &index, - &emptyIndex, - S2BooleanOperation::OpType::UNION, - options.booleanOperationOptions(), - options.layerOptions() - ); - - return List::create(Rcpp::XPtr(geography.release())); + std::unique_ptr geog_out = agg.Finalize(); + return List::create(Geography::MakeXPtr(std::move(geog_out))); } -// This approach to aggregation is slow but accurate. There is probably a more efficient way -// to accumulate geometries and/or re-use the layers vector but thus far I haven't figured -// out a way to make that work. // [[Rcpp::export]] List cpp_s2_union_agg(List geog, List s2options, bool naRm) { GeographyOperationOptions options(s2options); - GeographyOperationOptions::LayerOptions layerOptions = options.layerOptions(); - S2BooleanOperation::Options unionOptions = options.booleanOperationOptions(); - S2Builder::Options buillderOptions = options.builderOptions(); - - // using smart pointers here so that we can use swap() to - // use replace accumulatedIndex with index after each union - std::unique_ptr index = absl::make_unique(); - std::unique_ptr accumulatedIndex = absl::make_unique(); + s2geography::S2UnionAggregator agg(options.geographyOptions()); SEXP item; for (R_xlen_t i = 0; i < geog.size(); i++) { @@ -271,40 +91,17 @@ List cpp_s2_union_agg(List geog, List s2options, bool naRm) { if (item != R_NilValue) { Rcpp::XPtr feature(item); - - index->Clear(); - s2builderutil::LayerVector layers(3); - layers[0] = absl::make_unique(index.get(), layerOptions.pointLayerOptions); - layers[1] = absl::make_unique(index.get(), layerOptions.polylineLayerOptions); - layers[2] = absl::make_unique(index.get(), layerOptions.polygonLayerOptions); - - S2BooleanOperation booleanOp( - S2BooleanOperation::OpType::UNION, - s2builderutil::NormalizeClosedSet(std::move(layers)), - unionOptions - ); - - S2Error error; - if (!booleanOp.Build(*accumulatedIndex, *(feature->ShapeIndex()), &error)) { - stop(error.text()); - } - - accumulatedIndex.swap(index); + agg.Add(feature->Geog()); } } - std::unique_ptr geography = rebuildGeography( - accumulatedIndex.get(), - options.builderOptions(), - options.layerOptions() - ); - - return List::create(Rcpp::XPtr(geography.release())); + std::unique_ptr geog_out = agg.Finalize(); + return List::create(Geography::MakeXPtr(std::move(geog_out))); } // [[Rcpp::export]] List cpp_s2_centroid_agg(List geog, bool naRm) { - S2Point cumCentroid; + s2geography::CentroidAggregator agg; SEXP item; for (R_xlen_t i = 0; i < geog.size(); i++) { @@ -315,18 +112,17 @@ List cpp_s2_centroid_agg(List geog, bool naRm) { if (item != R_NilValue) { Rcpp::XPtr feature(item); - S2Point centroid = feature->Centroid(); - if (centroid.Norm2() > 0) { - cumCentroid += centroid.Normalize(); - } + agg.Add(feature->Geog()); } } + S2Point centroid = agg.Finalize(); + List output(1); - if (cumCentroid.Norm2() == 0) { - output[0] = Rcpp::XPtr(new PointGeography()); + if (centroid.Norm2() == 0) { + output[0] = Geography::MakeXPtr(Geography::MakePoint()); } else { - output[0] = Rcpp::XPtr(new PointGeography(cumCentroid.Normalize())); + output[0] = Geography::MakeXPtr(Geography::MakePoint(centroid)); } return output; @@ -336,7 +132,9 @@ List cpp_s2_centroid_agg(List geog, bool naRm) { List cpp_s2_rebuild_agg(List geog, List s2options, bool naRm) { GeographyOperationOptions options(s2options); - MutableS2ShapeIndex index; + s2geography::RebuildAggregator agg(options.geographyOptions()); + std::vector> geographies; + SEXP item; for (R_xlen_t i = 0; i < geog.size(); i++) { item = geog[i]; @@ -346,57 +144,12 @@ List cpp_s2_rebuild_agg(List geog, List s2options, bool naRm) { if (item != R_NilValue) { Rcpp::XPtr feature(item); - feature->BuildShapeIndex(&index); + agg.Add(feature->Geog()); } } - std::unique_ptr geography = rebuildGeography( - &index, - options.builderOptions(), - options.layerOptions() - ); - - return List::create(Rcpp::XPtr(geography.release())); -} - -std::vector findClosestPoints(S2ShapeIndex* index1, S2ShapeIndex* index2) { - // see http://s2geometry.io/devguide/s2closestedgequery.html section on Modeling Accuracy: - - // Find the edge from index2 that is closest to index1 - S2ClosestEdgeQuery query1(index1); - query1.mutable_options()->set_include_interiors(false); - S2ClosestEdgeQuery::ShapeIndexTarget target1(index2); - auto result1 = query1.FindClosestEdge(&target1); - - if (result1.edge_id() == -1) { - return std::vector(); - } - - // Get the edge from index1 (edge1) that is closest to index2. - S2Shape::Edge edge1 = query1.GetEdge(result1); - - // Now find the edge from index2 (edge2) that is closest to edge1. - S2ClosestEdgeQuery query2(index2); - query2.mutable_options()->set_include_interiors(false); - S2ClosestEdgeQuery::EdgeTarget target2(edge1.v0, edge1.v1); - auto result2 = query2.FindClosestEdge(&target2); - - // what if result2 has no edges? - if (result2.is_interior()) { - stop("S2ClosestEdgeQuery result is interior!"); - } - S2Shape::Edge edge2 = query2.GetEdge(result2); - - // Find the closest point pair on edge1 and edge2. - std::pair closest = S2::GetEdgePairClosestPoints( - edge1.v0, edge1.v1, - edge2.v0, edge2.v1 - ); - - std::vector pts(2); - pts[0] = closest.first; - pts[1] = closest.second; - return pts; + auto geog_out = agg.Finalize(); + return List::create(Geography::MakeXPtr(std::move(geog_out))); } // [[Rcpp::export]] @@ -404,12 +157,11 @@ List cpp_s2_closest_point(List geog1, List geog2) { class Op: public BinaryGeographyOperator { SEXP processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - std::vector pts = findClosestPoints(feature1->ShapeIndex(), feature2->ShapeIndex()); - - if (pts.size() == 0) { - return XPtr(new PointGeography()); + S2Point pt = s2geography::s2_closest_point(feature1->Index(), feature2->Index()); + if (pt.Norm2() == 0) { + return Geography::MakeXPtr(Geography::MakePoint()); } else { - return XPtr(new PointGeography(pts[0])); + return Geography::MakeXPtr(Geography::MakePoint(pt)); } } }; @@ -423,18 +175,28 @@ List cpp_s2_minimum_clearance_line_between(List geog1, List geog2) { class Op: public BinaryGeographyOperator { SEXP processFeature(XPtr feature1, XPtr feature2, R_xlen_t i) { - std::vector pts = findClosestPoints(feature1->ShapeIndex(), feature2->ShapeIndex()); + std::pair pts = s2geography::s2_minimum_clearance_line_between( + feature1->Index(), + feature2->Index() + ); + + if (pts.first.Norm2() == 0) { + return Geography::MakeXPtr(Geography::MakePoint()); + } + + std::vector vertices(2); + vertices[0] = pts.first; + vertices[1] = pts.second; - if (pts.size() == 0) { - return XPtr(new PolylineGeography()); - } else if (pts[0] == pts[1]) { - return XPtr(new PointGeography(pts)); + if (pts.first == pts.second) { + return Geography::MakeXPtr(Geography::MakePoint(std::move(vertices))); } else { + std::vector vertices(2); + vertices[0] = pts.first; + vertices[1] = pts.second; std::unique_ptr polyline = absl::make_unique(); - polyline->Init(pts); - std::vector> polylines(1); - polylines[0] = std::move(polyline); - return XPtr(new PolylineGeography(std::move(polylines))); + polyline->Init(vertices); + return Geography::MakeXPtr(Geography::MakePolyline(std::move(polyline))); } } }; @@ -447,11 +209,11 @@ List cpp_s2_minimum_clearance_line_between(List geog1, List geog2) { List cpp_s2_centroid(List geog) { class Op: public UnaryGeographyOperator { SEXP processFeature(XPtr feature, R_xlen_t i) { - S2Point centroid = feature->Centroid(); + S2Point centroid = s2geography::s2_centroid(feature->Geog()); if (centroid.Norm2() == 0) { - return XPtr(new PointGeography()); + return Geography::MakeXPtr(Geography::MakePoint()); } else { - return XPtr(new PointGeography(centroid.Normalize())); + return Geography::MakeXPtr(Geography::MakePoint(centroid.Normalize())); } } }; @@ -467,48 +229,11 @@ List cpp_s2_point_on_surface(List geog) { S2RegionCoverer coverer; SEXP processFeature(XPtr feature, R_xlen_t i) { - if (feature->IsEmpty()) { - return XPtr(new PointGeography()); - } - else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYGON) { - // Create Interior Covering of Polygon - S2CellUnion cellUnion; - cellUnion = coverer.GetInteriorCovering(*feature->Polygon()); - - // Take center of cell with smallest level (biggest) - int min_level = 31; - S2Point pt; - for (const S2CellId& id : cellUnion) { - if (id.level() < min_level) { - // Already normalized - // https://github.com/r-spatial/s2/blob/a323a74774d0526f1165e334a3d76a9da38316b9/src/s2/s2cell_id.h#L128 - pt = id.ToPoint(); - min_level = id.level(); - } - } - - return XPtr(new PointGeography(pt)); - } - // For point, return point closest to centroid - else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POINT) { - S2Point centroid = feature->Centroid(); - const std::vector* pts = feature->Point(); - - S1Angle nearest_dist = S1Angle::Infinity(); - S1Angle dist; - S2Point closest_pt; - for (const S2Point& pt : *pts) { - dist = S1Angle(pt, centroid); - if (dist < nearest_dist) { - nearest_dist = dist; - closest_pt = pt; - } - } - - return XPtr(new PointGeography(closest_pt)); - } - else { - stop("POLYLINE/COLLECTION type not supported in s2_point_on_surface()"); + S2Point result = s2geography::s2_point_on_surface(feature->Geog(), coverer); + if (result.Norm2() == 0) { + return Geography::MakeXPtr(Geography::MakePoint()); + } else { + return Geography::MakeXPtr(Geography::MakePoint(result)); } } }; @@ -521,8 +246,8 @@ List cpp_s2_point_on_surface(List geog) { List cpp_s2_boundary(List geog) { class Op: public UnaryGeographyOperator { SEXP processFeature(XPtr feature, R_xlen_t i) { - std::unique_ptr ptr = feature->Boundary(); - return XPtr(ptr.release()); + std::unique_ptr result = s2geography::s2_boundary(feature->Geog()); + return Geography::MakeXPtr(std::move(result)); } }; @@ -536,22 +261,20 @@ List cpp_s2_rebuild(List geog, List s2options) { public: Op(List s2options) { GeographyOperationOptions options(s2options); - this->options = options.builderOptions(); - this->layerOptions = options.layerOptions(); + this->options = options.geographyOptions(); } SEXP processFeature(XPtr feature, R_xlen_t i) { - std::unique_ptr ptr = rebuildGeography( - feature->ShapeIndex(), - this->options, - this->layerOptions + std::unique_ptr ptr = s2geography::s2_rebuild( + feature->Geog(), + this->options ); - return XPtr(ptr.release()); + + return Geography::MakeXPtr(std::move(ptr)); } private: - S2Builder::Options options; - GeographyOperationOptions::LayerOptions layerOptions; + s2geography::GlobalOptions options; }; Op op(s2options); @@ -564,102 +287,19 @@ List cpp_s2_unary_union(List geog, List s2options) { public: Op(List s2options) { GeographyOperationOptions options(s2options); - this->options = options.booleanOperationOptions(); - this->layerOptions = options.layerOptions(); + this->geographyOptions = options.geographyOptions(); } SEXP processFeature(XPtr feature, R_xlen_t i) { - // complex union only needed when a polygon is involved - bool simpleUnionOK = feature->IsEmpty() || - (feature->Dimension() < 2); - - // valid polygons that are not part of a collection can also use a - // simple union (common) - if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYGON) { - S2Error validationError; - if(!(feature->Polygon()->FindValidationError(&validationError))) { - simpleUnionOK = true; - } - } - - if (simpleUnionOK) { - MutableS2ShapeIndex emptyIndex; - - std::unique_ptr ptr = doBooleanOperation( - feature->ShapeIndex(), - &emptyIndex, - S2BooleanOperation::OpType::UNION, - this->options, - this->layerOptions - ); - - return XPtr(ptr.release()); - } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYGON) { - // If we've made it here we have an invalid polygon on our hands. A geography with - // invalid loops won't work with the S2BooleanOperation we will use to accumulate - // (i.e., union) valid polygons, so we need to rebuild each loop as its own polygon, - // splitting crossed edges along the way. - const S2Polygon* originalPoly = feature->Polygon(); - - // Not exposing these options as an argument (except snap function) - // because a particular combiation of them is required for this to work - S2Builder::Options builderOptions; - builderOptions.set_split_crossing_edges(true); - builderOptions.set_snap_function(this->options.snap_function()); - s2builderutil::S2PolygonLayer::Options layerOptions; - layerOptions.set_edge_type(S2Builder::EdgeType::UNDIRECTED); - layerOptions.set_validate(false); - - // Rebuild all loops as polygons using the S2Builder() - std::vector> loops; - for (int i = 0; i < originalPoly->num_loops(); i++) { - std::unique_ptr loop = absl::make_unique(); - S2Builder builder(builderOptions); - builder.StartLayer(absl::make_unique(loop.get())); - builder.AddShape(S2Loop::Shape(originalPoly->loop(i))); - S2Error error; - if (!builder.Build(&error)) { - throw GeographyOperatorException(error.text()); - } - - // Check if the builder created a polygon whose boundary contained more than - // half the earth (and invert it if so) - if (loop->GetArea() > (2 * M_PI)) { - loop->Invert(); - } - - loops.push_back(std::move(loop)); - } - - // Accumulate the union of outer loops (but difference of inner loops) - std::unique_ptr accumulatedPolygon = absl::make_unique(); - for (int i = 0; i < originalPoly->num_loops(); i++) { - std::unique_ptr polygonResult = absl::make_unique(); - - // Use original nesting to suggest if this loop should be unioned or diffed. - // For valid polygons loops are arranged such that the biggest loop is on the outside - // followed by holes such that the below strategy should work (since we are - // just iterating along the original loop structure) - if ((originalPoly->loop(i)->depth() % 2) == 0) { - polygonResult->InitToUnion(accumulatedPolygon.get(), loops[i].get()); - } else { - polygonResult->InitToDifference(accumulatedPolygon.get(), loops[i].get()); - } - - accumulatedPolygon.swap(polygonResult); - } - - return XPtr(new PolygonGeography(std::move(accumulatedPolygon))); - } else { - // This is a less common case (mixed dimension output that includes a polygon). - // In the absence of a clean solution, saving this battle for another day. - throw GeographyOperatorException("Unary union for collections is not implemented"); - } + std::unique_ptr geog_out = + s2geography::s2_unary_union(feature->Index(), this->geographyOptions); + return Geography::MakeXPtr(std::move(geog_out)); } private: S2BooleanOperation::Options options; GeographyOperationOptions::LayerOptions layerOptions; + s2geography::GlobalOptions geographyOptions; }; Op op(s2options); @@ -677,19 +317,22 @@ List cpp_s2_interpolate_normalized(List geog, NumericVector distanceNormalized) return R_NilValue; } - if (feature->IsCollection()) { - throw GeographyOperatorException("`x` must be a simple geography"); + if (s2geography::s2_is_empty(feature->Geog())) { + return Geography::MakeXPtr(Geography::MakePoint()); } - if (feature->IsEmpty()) { - return R_NilValue; + if (s2geography::s2_is_collection(feature->Geog())) { + throw GeographyOperatorException("`x` must be a simple geography"); + } else if (feature->Geog().dimension() != 1) { + throw GeographyOperatorException("`x` must be a polyline"); } - if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYLINE) { - S2Point point = feature->Polyline()->at(0)->Interpolate(this->distanceNormalized[i]); - return XPtr(new PointGeography(point)); + S2Point point = s2geography::s2_interpolate_normalized(feature->Geog(), this->distanceNormalized[i]); + + if (point.Norm2() == 0) { + return Geography::MakeXPtr(Geography::MakePoint()); } else { - throw GeographyOperatorException("`x` must be a polyline geography"); + return Geography::MakeXPtr(Geography::MakePoint(point)); } } }; @@ -714,7 +357,7 @@ List cpp_s2_buffer_cells(List geog, NumericVector distance, int maxCells, int mi SEXP processFeature(XPtr feature, R_xlen_t i) { S2ShapeIndexBufferedRegion region; - region.Init(feature->ShapeIndex(), S1ChordAngle::Radians(this->distance[i])); + region.Init(&feature->Index().ShapeIndex(), S1ChordAngle::Radians(this->distance[i])); S2CellUnion cellUnion; cellUnion = coverer.GetCovering(region); @@ -722,7 +365,7 @@ List cpp_s2_buffer_cells(List geog, NumericVector distance, int maxCells, int mi std::unique_ptr polygon = absl::make_unique(); polygon->InitToCellUnionBorder(cellUnion); - return XPtr(new PolygonGeography(std::move(polygon))); + return Geography::MakeXPtr(Geography::MakePolygon(std::move(polygon))); } }; @@ -730,46 +373,13 @@ List cpp_s2_buffer_cells(List geog, NumericVector distance, int maxCells, int mi return op.processVector(geog); } - -class ConvexHullGeographyQuery: public S2ConvexHullQuery { -public: - - void AddGeography(Geography* feature) { - if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYGON) { - const S2Polygon* pol = feature->Polygon(); - this->AddPolygon(*pol); - } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POINT) { - const std::vector* pts = feature->Point(); - for(const auto& pt: *pts) { - this->AddPoint(pt); - } - } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_POLYLINE) { - const std::vector>* lins = feature->Polyline(); - for(const auto& lin: *lins) { - this->AddPolyline(*lin); - } - } else if (feature->GeographyType() == Geography::Type::GEOGRAPHY_COLLECTION) { - const std::vector>* features = feature->CollectionFeatures(); - for (const auto& feat: *features) { - this->AddGeography(feat.get()); - } - } - } - - std::unique_ptr GetConvexHullPolygon() { - return absl::make_unique(this->GetConvexHull()); - } -}; - // [[Rcpp::export]] List cpp_s2_convex_hull(List geog) { class Op: public UnaryGeographyOperator { SEXP processFeature(XPtr feature, R_xlen_t i) { - ConvexHullGeographyQuery convexHullQuery; - convexHullQuery.AddGeography(feature); - - std::unique_ptr outP = convexHullQuery.GetConvexHullPolygon(); - return XPtr(new PolygonGeography(std::move(outP))); + std::unique_ptr geog_out = + s2geography::s2_convex_hull(feature->Geog()); + return Geography::MakeXPtr(std::move(geog_out)); } }; @@ -780,8 +390,8 @@ List cpp_s2_convex_hull(List geog) { // [[Rcpp::export]] List cpp_s2_convex_hull_agg(List geog, bool naRm) { - // create the convex hull query - ConvexHullGeographyQuery convexHullQuery; + s2geography::S2ConvexHullAggregator agg; + SEXP item; for (R_xlen_t i = 0; i < geog.size(); i++) { item = geog[i]; @@ -791,11 +401,9 @@ List cpp_s2_convex_hull_agg(List geog, bool naRm) { if (item != R_NilValue) { XPtr feature(item); - convexHullQuery.AddGeography(feature); + agg.Add(feature->Geog()); } } - std::unique_ptr outP = convexHullQuery.GetConvexHullPolygon(); - XPtr outG(new PolygonGeography(std::move(outP))); - return List::create(outG); + return List::create(Geography::MakeXPtr(agg.Finalize())); } diff --git a/src/s2-c-api.cpp b/src/tessellator.cpp similarity index 98% rename from src/s2-c-api.cpp rename to src/tessellator.cpp index 28d55e89..3feb4445 100644 --- a/src/s2-c-api.cpp +++ b/src/tessellator.cpp @@ -4,13 +4,13 @@ #include "s2/s2projections.h" #include "s2/s2edge_tessellator.h" +// capi typedef start typedef struct s2_projection_t s2_projection_t; typedef struct s2_tessellator_t s2_tessellator_t; +// capi typedef end -#ifdef __cplusplus extern "C" { -#endif - +// capi func start s2_projection_t* s2_projection_create_plate_carree(double scale); s2_projection_t* s2_projection_create_mercator(double max_x); void s2_projection_destroy(s2_projection_t* projection); @@ -26,10 +26,8 @@ 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 +// capi func end } -#endif s2_projection_t* s2_projection_create_plate_carree(double scale) { return (s2_projection_t*) new S2::PlateCarreeProjection(scale); @@ -69,7 +67,7 @@ class TessellatorWrapper { TessellatorWrapper(s2_projection_t* projection, double tolerance_radians): tessellator((S2::Projection*) projection, S1Angle::Radians(tolerance_radians)), has_r2_last(false), has_s2_last(false) {} - + void reset() { s2points.clear(); r2points.clear(); diff --git a/src/tessellator.h b/src/tessellator.h new file mode 100644 index 00000000..3e06c161 --- /dev/null +++ b/src/tessellator.h @@ -0,0 +1,27 @@ + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct s2_projection_t s2_projection_t; +typedef struct s2_tessellator_t s2_tessellator_t; + +s2_projection_t* s2_projection_create_plate_carree(double scale); +s2_projection_t* s2_projection_create_mercator(double max_x); +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 diff --git a/src/wk-c-utils.c b/src/wk-c-utils.c index 4c24ca5e..8e6f27ec 100644 --- a/src/wk-c-utils.c +++ b/src/wk-c-utils.c @@ -3,33 +3,7 @@ #include #include #include "wk-v1.h" - -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); -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 +#include "tessellator.h" // Expose projections as external pointers so that they can theoretically be generated // by other packages. @@ -348,7 +322,7 @@ SEXP c_s2_coord_filter_new(SEXP handler_xptr, SEXP projection_xptr, SEXP unproje } else { coord_filter->tessellator = NULL; } - + coord_filter->unproject = LOGICAL(unproject)[0]; if (coord_filter->unproject) { handler->coord = &s2_coord_filter_coord_unproject; diff --git a/src/wk-geography.h b/src/wk-geography.h deleted file mode 100644 index 1ca75ad4..00000000 --- a/src/wk-geography.h +++ /dev/null @@ -1,168 +0,0 @@ - -#ifndef WK_GEOGRAPHY_H -#define WK_GEOGRAPHY_H - -#include - -#include "wk/rcpp-io.hpp" -#include "wk/reader.hpp" -#include "wk/geometry-handler.hpp" - -#include "point-geography.h" -#include "polyline-geography.h" -#include "polygon-geography.h" -#include "geography-collection.h" - -#include -#include "geography.h" - -#define CODE_HAS_BUILD_ERROR 3938829 - - -class WKGeographyWriter: public WKGeometryHandler { -public: - Rcpp::List output; - R_xlen_t featureId; - - Rcpp::IntegerVector problemId; - Rcpp::CharacterVector problems; - - WKGeographyWriter(R_xlen_t size): - output(size), - builder(nullptr), - oriented(false), - check(true) {} - - void setOriented(bool oriented) { - this->oriented = oriented; - } - - void setCheck(bool check) { - this->check = check; - } - - void nextFeatureStart(size_t featureId) { - this->builder = std::unique_ptr(nullptr); - this->featureId = featureId; - } - - void nextNull(size_t featureId) { - this->output[featureId] = R_NilValue; - } - - void nextGeometryStart(const WKGeometryMeta& meta, uint32_t partId) { - if (!this->builder) { - switch (meta.geometryType) { - case WKGeometryType::Point: - case WKGeometryType::MultiPoint: - this->builder = absl::make_unique(); - break; - case WKGeometryType::LineString: - case WKGeometryType::MultiLineString: - this->builder = absl::make_unique(); - break; - case WKGeometryType::Polygon: - case WKGeometryType::MultiPolygon: - this->builder = absl::make_unique( - this->oriented, - this->check - ); - break; - case WKGeometryType::GeometryCollection: - this->builder = absl::make_unique( - this->oriented, - this->check - ); - break; - default: - std::stringstream err; - err << "Unknown geometry type in geography builder: " << meta.geometryType; - this->addProblem(err.str()); - throw WKParseException(CODE_HAS_BUILD_ERROR); - } - } - - this->builder->nextGeometryStart(meta, partId); - } - - void nextLinearRingStart(const WKGeometryMeta& meta, uint32_t size, uint32_t ringId) { - this->builder->nextLinearRingStart(meta, size, ringId); - } - - void nextCoordinate(const WKGeometryMeta& meta, const WKCoord& coord, uint32_t coordId) { - this->builder->nextCoordinate(meta, coord, coordId); - } - - void nextLinearRingEnd(const WKGeometryMeta& meta, uint32_t size, uint32_t ringId) { - try { - this->builder->nextLinearRingEnd(meta, size, ringId); - } catch (WKParseException& e) { - this->addProblem(e.what()); - throw WKParseException(CODE_HAS_BUILD_ERROR); - } - } - - void nextGeometryEnd(const WKGeometryMeta& meta, uint32_t partId) { - this->builder->nextGeometryEnd(meta, partId); - } - - void nextFeatureEnd(size_t featureId) { - if (this->builder) { - try { - std::unique_ptr feature = builder->build(); - this->output[featureId] = Rcpp::XPtr(feature.release()); - } catch (WKParseException& e) { - this->addProblem(e.what()); - throw WKParseException(CODE_HAS_BUILD_ERROR); - } - } - } - - bool nextError(WKParseException& error, size_t featureId) { - if (error.code() == CODE_HAS_BUILD_ERROR) { - this->output[featureId] = R_NilValue; - return true; - } else { - return false; - } - - this->nextFeatureEnd(featureId); - return true; - } - -private: - std::unique_ptr builder; - bool oriented; - bool check; - - void addProblem(std::string what) { - problemId.push_back(this->featureId); - problems.push_back(what); - } -}; - - -class WKGeographyReader: public WKReader { -public: - - WKGeographyReader(WKRcppSEXPProvider& provider): - WKReader(provider), provider(provider) {} - - void readFeature(size_t featureId) { - this->handler->nextFeatureStart(featureId); - - if (this->provider.featureIsNull()) { - this->handler->nextNull(featureId); - } else { - Rcpp::XPtr geography(this->provider.feature()); - geography->Export(handler, WKReader::PART_ID_NONE); - } - - this->handler->nextFeatureEnd(featureId); - } - -private: - WKRcppSEXPProvider& provider; -}; - -#endif diff --git a/tests/testthat/test-s2-accessors.R b/tests/testthat/test-s2-accessors.R index e929f1d2..cac91f00 100644 --- a/tests/testthat/test-s2-accessors.R +++ b/tests/testthat/test-s2-accessors.R @@ -29,7 +29,7 @@ test_that("s2_is_valid() works", { # invalid "LINESTRING (0 0, 0 0, 1 1)", "POLYGON ((0 0, 0 1, 1 0, 0 0, 0 0))", - "GEOMETRYCOLLECTION (POLYGON ((0 0, 0 1, 1 0, 0 0, 0 0)))", + "GEOMETRYCOLLECTION (POINT (0 1), POLYGON ((0 0, 0 1, 1 0, 0 0, 0 0)))", NA ) ), @@ -50,7 +50,7 @@ test_that("s2_is_valid_detail() works", { # invalid "LINESTRING (0 0, 0 0, 1 1)", "POLYGON ((0 0, 0 1, 1 0, 0 0, 0 0))", - "GEOMETRYCOLLECTION (POLYGON ((0 0, 0 1, 1 0, 0 0, 0 0)))", + "GEOMETRYCOLLECTION (POINT (0 1), POLYGON ((0 0, 0 1, 1 0, 0 0, 0 0)))", NA ) ), @@ -81,6 +81,10 @@ test_that("s2_num_points works", { expect_identical(s2_num_points("POINT EMPTY"), 0L) expect_identical(s2_num_points("LINESTRING (0 0, 1 1)"), 2L) expect_identical(s2_num_points("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))"), 4L) + expect_identical( + s2_num_points("GEOMETRYCOLLECTION (POINT (0 1), LINESTRING (0 0, 1 1))"), + 3L + ) }) test_that("s2_is_empty works", { @@ -104,6 +108,11 @@ test_that("s2_area works", { expect_true( abs(s2_area("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))", radius = 180 / pi) - 100) < 0.27 ) + + expect_identical( + s2_area("POLYGON ((0 0, 90 0, 0 90, 0 0))"), + s2_area("GEOMETRYCOLLECTION(POLYGON ((0 0, 90 0, 0 90, 0 0)))") + ) }) test_that("s2_length works", { @@ -134,8 +143,8 @@ test_that("s2_x and s2_y works", { expect_identical(s2_y(NA_character_), NA_real_) expect_equal(s2_x("POINT (-64 45)"), -64) expect_equal(s2_y("POINT (-64 45)"), 45) - expect_identical(s2_x("POINT EMPTY"), NA_real_) - expect_identical(s2_y("POINT EMPTY"), NA_real_) + expect_identical(s2_x("POINT EMPTY"), NaN) + expect_identical(s2_y("POINT EMPTY"), NaN) expect_error(s2_x("LINESTRING EMPTY"), "Can't compute") expect_error(s2_y("LINESTRING EMPTY"), "Can't compute") expect_error(s2_x("POLYGON EMPTY"), "Can't compute") @@ -160,17 +169,17 @@ test_that("s2_project() and s2_project_normalized() work", { c(0, 0.25, 0.75, 1, NA_real_, NA_real_) ) - expect_error( + expect_identical( s2_project_normalized("POINT (0 1)", "POINT (0 1)"), - "must be a polyline" + NaN ) - expect_error( + expect_identical( s2_project_normalized("LINESTRING (0 1, 1 1)", "LINESTRING (0 1, 1 1)"), - "must be a point" + NaN ) - expect_error( + expect_identical( s2_project_normalized("LINESTRING (0 1, 1 1)", "MULTIPOINT (0 1, 1 1)"), - "must both be simple geographies" + NaN ) }) diff --git a/tests/testthat/test-s2-geography.R b/tests/testthat/test-s2-geography.R index 106a672e..80d86af5 100644 --- a/tests/testthat/test-s2-geography.R +++ b/tests/testthat/test-s2-geography.R @@ -38,7 +38,7 @@ test_that("s2_geography vectors can be created from s2_lnglat and s2_point", { test_that("s2_geography vectors can be created from WKB and WKT", { wkb_point <- wk::as_wkb(wk::wkt("POINT (-64 45)", geodesic = TRUE)) - expect_output(print(as_s2_geography(wkb_point)), "") + expect_output(print(as_s2_geography(wkb_point)), "POINT \\(-64 45\\)") expect_error( as_s2_geography(wk::as_wkb("LINESTRING (0 0, 1 1)")), "Cartesian wkb\\(\\)" @@ -51,7 +51,7 @@ test_that("s2_geography vectors can be created from WKB and WKT", { expect_silent(as_s2_geography(wk::as_wkb("MULTIPOINT (0 1)"))) wkt_point <- wk::as_wkt(wk::wkt("POINT (-64 45)", geodesic = TRUE)) - expect_output(print(as_s2_geography(wkt_point)), "") + expect_output(print(as_s2_geography(wkt_point)), "POINT \\(-64 45\\)") expect_error( as_s2_geography(wk::wkt("LINESTRING (0 0, 1 1)")), "Cartesian wkt\\(\\)" @@ -64,8 +64,8 @@ test_that("s2_geography vectors can be created from WKB and WKT", { expect_silent(as_s2_geography(wk::wkt("MULTIPOINT (0 1)"))) # also test other classes commonly used to signify WKB or WKT - expect_output(print(as_s2_geography(structure(wkb_point, class = "WKB")), "")) - expect_output(print(as_s2_geography(structure(wkb_point, class = "blob")), "")) + expect_output(print(as_s2_geography(structure(wkb_point, class = "WKB")), "POINT \\(-64 45\\)")) + expect_output(print(as_s2_geography(structure(wkb_point, class = "blob")), "POINT \\(-64 45\\)")) }) test_that("s2_geography can be exported to WKB/WKT", { @@ -82,34 +82,34 @@ test_that("s2_geography can be exported to WKB/WKT", { }) test_that("s2_geography vectors can be created from wkt", { - expect_output(print(as_s2_geography("POINT (-64 45)")), "") - expect_output(print(as_s2_geography("POINT EMPTY")), "") + expect_output(print(as_s2_geography("POINT (-64 45)")), "POINT \\(-64 45\\)") + expect_output(print(as_s2_geography("POINT EMPTY")), "POINT EMPTY") expect_output( print(as_s2_geography("MULTIPOINT ((-64 45), (30 10))")), - "" + "MULTIPOINT \\(\\(-64 45\\), \\(30 10\\)\\)" ) expect_output( print(as_s2_geography("LINESTRING (-64 45, 0 0)")), - "" + "LINESTRING \\(-64 45, 0 0\\)" ) expect_output( print(as_s2_geography("LINESTRING EMPTY")), - "" + "LINESTRING EMPTY" ) expect_output( print(as_s2_geography("MULTILINESTRING ((-64 45, 0 0), (0 1, 2 3))")), - "" + "MULTILINESTRING \\(\\(-64 45, 0 0), \\(0 1, 2 3\\)\\)" ) - expect_output(print(as_s2_geography("POLYGON EMPTY"), "")) + expect_output(print(as_s2_geography("POLYGON EMPTY"), "POLYGON EMPTY")) expect_output( print(as_s2_geography("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))")), - "" + "GEOMETRYCOLLECTION \\(POINT \\(-64 45\\)\\)" ) expect_output( @@ -147,7 +147,7 @@ test_that("s2_geography vectors can be created from wkt", { ) ) - expect_output(print(as_s2_geography("GEOMETRYCOLLECTION EMPTY")), "") + expect_output(print(as_s2_geography("GEOMETRYCOLLECTION EMPTY")), "GEOMETRYCOLLECTION EMPTY") }) test_that("empty points are empty when imported from WKB", { @@ -168,7 +168,7 @@ test_that("nested ring depths are correctly exported", { )"), max_coords = 100 ), - "\\(20 35, 10 30, 10 10, 30 5, 45 20, 20 35\\), \\(30 20, 20 15, 20 25, 30 20" + "\\(20 35, 10 30, 10 10, 30 5, 45 20, 20 35\\), \\(30 20, 20 15, 20 25" ) # polygon with a hole in a hole! @@ -178,13 +178,13 @@ test_that("nested ring depths are correctly exported", { ((40 40, 20 45, 45 30, 40 40)), ( (20 35, 10 30, 10 10, 30 5, 45 20, 20 35), - (30 20, 20 15, 20 25, 30 20), - (27 21, 21 21, 21 16, 27 21) - ) + (30 20, 20 15, 20 25, 30 20) + ), + ((27 21, 21 21, 21 16, 27 21)) )"), max_coords = 100 ), - "30 20, 20 15, 20 25, 30 20\\), \\(27 21, 21 21, 21 16, 27 21" + "\\(\\(27 21, 21 21, 21 16, 27 21\\)\\)\\)" ) }) diff --git a/tests/testthat/test-s2-transformers.R b/tests/testthat/test-s2-transformers.R index 2fa99f45..25a1ae0d 100644 --- a/tests/testthat/test-s2-transformers.R +++ b/tests/testthat/test-s2-transformers.R @@ -279,7 +279,7 @@ test_that("s2_union(x) errors for the case of mixed dimension collections", { s2_union( c("GEOMETRYCOLLECTION(POLYGON ((-10 -10, -10 10, 10 10, 10 -10, -10 -10)), LINESTRING (0 -20, 0 20))") ), - "Unary union for collections is not implemented" + "for multidimensional collections not implemented" ) }) @@ -388,17 +388,36 @@ test_that("s2_union_agg() works", { # NULL handling expect_identical( - s2_coverage_union_agg(c("POINT (30 10)", NA), na.rm = FALSE), + s2_union_agg(c("POINT (30 10)", NA), na.rm = FALSE), as_s2_geography(NA_character_) ) expect_wkt_equal( - s2_coverage_union_agg(character()), + s2_union_agg(character()), as_s2_geography("GEOMETRYCOLLECTION EMPTY") ) expect_wkt_equal( - s2_coverage_union_agg(c("POINT (30 10)", NA), na.rm = TRUE), + s2_union_agg(c("POINT (30 10)", NA), na.rm = TRUE), "POINT (30 10)" ) + + # make sure this works on polygons since they are handled differently than + # points and linestrings + expect_equal( + s2_area(s2_union_agg(s2_data_countries())), + sum(s2_area(s2_union_agg(s2_data_countries()))) + ) + + # check non-polygons and polygons together + points_and_poly <- s2_union_agg( + c( + s2_data_countries(), + s2_data_cities() + ) + ) + + points <- s2_rebuild(points_and_poly, options = s2_options(dimensions = "point")) + poly <- s2_rebuild(points_and_poly, options = s2_options(dimensions = "polygon")) + expect_false(any(s2_intersects(points, poly))) }) test_that("s2_rebuild_agg() works", { diff --git a/tests/testthat/test-wk-utils.R b/tests/testthat/test-wk-utils.R index f8d657bb..a2d884e5 100644 --- a/tests/testthat/test-wk-utils.R +++ b/tests/testthat/test-wk-utils.R @@ -184,3 +184,34 @@ test_that("mercator projection works", { wk::xy(c(0, 20037508), 0) ) }) + +test_that("wk_handle() for s2_geography works", { + for (name in names(s2_data_example_wkt)) { + geog <- wk::wk_handle( + s2_data_example_wkt[[name]], + s2_geography_writer() + ) + + geog2 <- wk::wk_handle( + geog, + s2_geography_writer(check = TRUE, oriented = TRUE) + ) + + expect_equal(wk::wk_coords(geog), wk::wk_coords(geog2)) + } +}) + +test_that("the s2_geography_writer() works", { + # nc has some rings that get reordered by this operation + for (name in setdiff(names(s2_data_example_wkt), "nc")) { + geog <- wk::wk_handle( + s2_data_example_wkt[[name]], + s2_geography_writer() + ) + + expect_equal( + wk::wk_coords(as_wkt(geog))[c("x", "y")], + wk::wk_coords(s2_data_example_wkt[[name]])[c("x", "y")] + ) + } +})