geos_coord_seq.Rmd

```{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;
}
```