This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
```{r setup, include = FALSE} | |
library(sf) | |
library(wk) | |
library(cpp11) | |
Sys.setenv( | |
PKG_CXXFLAGS = paste0("-I", getwd()) | |
) | |
``` | |
```{r} | |
library(arrowvctrs) | |
geo_line_schema <- arrow_schema( | |
"+l", | |
children = list( | |
arrow_schema( | |
"+s", | |
children = list( | |
arrow_schema("g", name = "x"), | |
arrow_schema("g", name = "y") | |
) | |
) | |
) | |
) | |
geo_line <- arrow_vctr( | |
geo_line_schema, | |
arrow_array( | |
buffers = c(0L, 2L, 4L), | |
length = 2L, | |
children = list( | |
arrow_array( | |
length = 4L, | |
children = list( | |
arrow_array(c(0, 5, 0, 10), length = 4L), | |
arrow_array(c(0, 10, 5, 0), length = 4L) | |
) | |
) | |
) | |
) | |
) | |
``` | |
```{r, include=FALSE} | |
Sys.setenv( | |
PKG_CXXFLAGS = paste0("-I", getwd(), " -I/opt/homebrew/Cellar/geos/3.9.1/include"), | |
PKG_LIBS = '-L/opt/homebrew/Cellar/geos/3.9.1/lib -lgeos-3.9.1 -lgeos' | |
) | |
``` | |
```{cpp11} | |
#include "cpp11.hpp" | |
#include "abi.h" | |
#include <geos/geom/GeometryFactory.h> | |
#include <geos/geom/LineString.h> | |
#include <geos/geom/CoordinateFilter.h> | |
#include <geos/geom/CoordinateSequence.h> | |
#include <memory.h> | |
using namespace cpp11; | |
using namespace geos; | |
using namespace geos::geom; | |
class ArrowCoordinateSequence: public CoordinateSequence { | |
public: | |
ArrowCoordinateSequence(struct ArrowArray* geo_line, int size, int base_offset): | |
size(size), base_offset(base_offset) { | |
Rprintf("ArrowCoordinateSequence\n"); | |
this->x = (double*) geo_line->children[0]->children[0]->buffers[0]; | |
this->y = (double*) geo_line->children[0]->children[1]->buffers[0]; | |
this->toVector(this->coords); | |
} | |
virtual std::size_t getSize() const { | |
Rprintf("getSize\n"); | |
return this->size; | |
} | |
virtual const Coordinate& getAt(std::size_t i) const { | |
Rprintf("getAt %d\n", i); | |
return this->coords[i]; | |
return coord; | |
} | |
virtual void getAt(std::size_t i, Coordinate& c) const { | |
Rprintf("getAtWrite %d\n", i); | |
c.x = this->x[this->base_offset + i]; | |
c.y = this->y[this->base_offset + i]; | |
} | |
virtual double getOrdinate(std::size_t index, std::size_t ordinateIndex) const { | |
Rprintf("getOrdinate\n"); | |
if (ordinateIndex == 0) { | |
return this->x[this->base_offset + index]; | |
} else { | |
return this->y[this->base_offset + index]; | |
} | |
} | |
virtual void toVector(std::vector<Coordinate>& coords) const { | |
Rprintf("toVector\n"); | |
double x, y; | |
for (int i = 0; i < this->size; i++) { | |
x = this->x[this->base_offset + i]; | |
y = this->y[this->base_offset + i]; | |
coords.push_back(Coordinate(x, y)); | |
} | |
} | |
virtual bool isEmpty() const { | |
Rprintf("isEmpty\n"); | |
return this->size == 0; | |
} | |
virtual void apply_ro(CoordinateFilter* filter) const { | |
Rprintf("apply_ro\n"); | |
Coordinate c; | |
for (int i = 0; i < this->size; i++) { | |
c.x = this->x[this->base_offset + i]; | |
c.y = this->y[this->base_offset + i]; | |
filter->filter_ro(&c); | |
} | |
} | |
virtual std::size_t getDimension() const { | |
Rprintf("getDimension\n"); | |
return 2; | |
} | |
virtual void setPoints(const std::vector<Coordinate>& v) { | |
stop("not supported"); | |
} | |
virtual std::unique_ptr<CoordinateSequence> clone() const { | |
stop("not supported"); | |
} | |
virtual void apply_rw(const CoordinateFilter* filter) { | |
stop("not supported"); | |
} | |
virtual void setAt(const Coordinate& c, std::size_t pos) { | |
stop("not supported"); | |
} | |
virtual void setOrdinate(std::size_t index, std::size_t ordinateIndex, double value) { | |
stop("not supported"); | |
} | |
private: | |
int size; | |
int base_offset; | |
double* x; | |
double* y; | |
std::vector<Coordinate> coords; | |
// if you just return a reference to this coordinate in getAt() | |
// then intersects() gives you the wrong answer | |
mutable Coordinate coord; | |
}; | |
[[cpp11::register]] | |
logicals arrow_line_array_geos_intersects(sexp geo_line_array1, sexp geo_line_array2) { | |
struct ArrowArray* geo_line1 = (struct ArrowArray*) R_ExternalPtrAddr(geo_line_array1); | |
int* offsets1 = (int*) geo_line1->buffers[0]; | |
struct ArrowArray* geo_line2 = (struct ArrowArray*) R_ExternalPtrAddr(geo_line_array2); | |
int* offsets2 = (int*) geo_line2->buffers[0]; | |
if (geo_line1->length != geo_line2->length) { | |
stop("Both arrays must be of equal length"); | |
} | |
const GeometryFactory* factory = GeometryFactory::getDefaultInstance(); | |
writable::logicals result(geo_line1->length); | |
for (int i = 0; i < geo_line1->length; i++) { | |
int line_size1 = offsets1[i + 1] - offsets1[i]; | |
ArrowCoordinateSequence* cs1 = new ArrowCoordinateSequence(geo_line1, line_size1, offsets1[i]); | |
std::unique_ptr<LineString> ls1(factory->createLineString(cs1)); | |
int line_size2 = offsets2[i + 1] - offsets2[i]; | |
ArrowCoordinateSequence* cs2 = new ArrowCoordinateSequence(geo_line2, line_size2, offsets2[i]); | |
std::unique_ptr<LineString> ls2(factory->createLineString(cs2)); | |
result[i] = ls1->intersects(ls2.get()); | |
} | |
return result; | |
} | |
``` |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment