Skip to content

Instantly share code, notes, and snippets.

@paleolimbot
Last active Sep 22, 2021
Embed
What would you like to do?
```{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