SpatialPolygons to geojson with Rcpp

SpatialPolygons_to_geojson.R

# convert SpatialPolygons/SpatialPolygonsDataFrame to geoJSON
# 2016-06 @javrucebo
#
# libraries used: sp, Rcpp, jsonlite, stringi
#
# FUNCTION SpatialPolygons_to_geojson
# parameter:
# sp: SpatialPolygons or SpatialPolygonsDataFrame
# precision: precision of numerical variables (default: 15)

# the cpp function take as parameter:
# sp: SpatialPolygons
# data: an array of JSON string with the data per feature in sp
# use_id: a boolean indicating whether to use the ID's from sp
#         (in case they are all numeric)
# precision: the precision for numeric values, default (14)
Rcpp::cppFunction('std::string sp_to_geojson_rcpp(Rcpp::S4 sp,
                                   Rcpp::CharacterVector data,
                                   bool use_id,
                                   int precision = 15) {
  // output string
  std::string out;

  /**********************************************************
  sp has a structure like:
  sp@polygons = list(Ps1, Ps2, ...)
  Psx = class Polygons = Polygons(list(Pxa, Pxb, ...), ID=id)
  Pxy = class Polygon = Polygon(coords)
  **********************************************************/

  Rcpp::List SpatialPolygons = sp.slot("polygons");

  out += "{\\n\\"type\\": \\"FeatureCollection\\",\\n                                                                                \\n\\"features\\": [\\n";

  int nPolygons = SpatialPolygons.size();
  int nrow_data = data.size();
  if (nPolygons != nrow_data) {
    return(std::string("ERROR: polygons and data differ in length"));
  }

  // check if we use "Polygon" or "Multipolygon" in JSON string.
  // Heuristics show that rgdal uses Multipolygon if one of the Polygons consists
  // of multiple rings
  bool multi = FALSE;
  for (int i=0; i < nPolygons; i++) {
    Rcpp::S4 Polygons = SpatialPolygons[i];
    Rcpp::List poly_list = Polygons.slot("Polygons");
    if (poly_list.size()>1) { // we found one multi-ring, so we are done.
       multi = TRUE;
       break;
    }
  }

  // Loop over slot "Polygons"
  for (int k=0; k<nPolygons; k++) {
    Rcpp::S4 Polygons = SpatialPolygons[k];
    Rcpp::List poly_list = Polygons.slot("Polygons");

    // if we do not use original (numeric) id, then take a sequence from loop
    std::string id;
    if (use_id) {
      Rcpp::String sid = Polygons.slot("ID");
      id = std::string(sid);
    } else {
      std::ostringstream s;
      s << k;
      id = s.str();
    }
    out += "{ \\"type\\": \\"Feature\\", \\"id\\": ";
    out += id;
    out += ", \\"properties\\": ";
    // write the data part
    out += data[k];
    out += ", \\"geometry\\": { \\"type\\": \\"";
    // MultiPolygon?
    if (multi) {
      out += "Multi";
    }
    out += "Polygon\\", \\"coordinates\\": [ ";
    if (multi) { // for multi add additional level or brackets
      out += "[ ";
    }

    // Write out all Polygons which are part of the feature
    int nPolygon = poly_list.size();
    for(int i=0; i<nPolygon; i++) {
      std::ostringstream s;
      Rcpp::S4 Polygon = poly_list[i];
      SEXP coords = Polygon.slot("coords");
      Rcpp::NumericMatrix x(coords);
      int nrow = x.nrow(), ncol = x.ncol();
      // if next ring is a hole we need to delete last closing/opening
      // square brackets and put the hole right after last ring.
      bool hole = Polygon.slot("hole");
      if (hole) {
        out.resize (out.size () - 6);
        out += ", ";
      }
      // starting bracket;
      out += "[ ";

      // stringstream to hold the coordinates
      s << std::setprecision(precision);
      // first to second last row
      for (int j = 0; j < nrow-1; j++) {
        s << "[ " << x(j,0) << ", " << x(j,1) << " ], ";
      }
      // last row without trailing comma
      s << "[ " << x(nrow-1,0) << ", " << x(nrow-1,1) << " ]";

      out += s.str();

      // closing bracket
      out += " ]";

      if (i < nPolygon -1) { // only if not last Polygon
        out += " ], [ ";
      }
    }

    if (multi) { // for multi close brackets
      out += " ]";
    }
    out += " ] } }";

    if (k < nPolygons-1) { // only inf not last Polygons
       out += ",\\n";
    }
  }

  out += "\\n]\\n}\\n";
  return out;
}')

SpatialPolygons_to_geojson <- function(sp, precision=15) {
  # if 'sp' is SpatialPolygonsDataFrame, we construct the JSON part of the data
  # here. otherwise dummy data
  if (inherits(sp, "SpatialPolygonsDataFrame")) {
    data <- jsonlite::toJSON(sp@data, na='null')
    data <- stringi::stri_sub(data,2,stri_length(data)-1)
    data <- stringi::stri_replace_all_regex(data,"\\},\\{","}@@@{")
    data <- stringi::stri_split_regex(data,"@@@")[[1]]
  } else if (inherits(sp, "SpatialPolygons")) {
    data <- rep('{ "dummy": 0.0 }', length(sp))
  } else {
    warning("Parameter 'sp' needs to be from classes SpatialPolygons(DataFrame)")
    return(NULL)
  }
  # currently we disregard the proj4string, as for many projections rgdal also
  # diregards informaiton.
  # TODO: transfer projection strings also to geoJSON
  # see http://geojson.org/geojson-spec.html#coordinate-reference-system-objects
  sp::proj4string(sp) <- ""
  # check whether only numeric id's are used for SpatialPolygons.
  # if not, for geojson we will create a new sequence of numerical id's
  use_ID <- !anyNA(as.numeric(sapply(slot(sp, "polygons"), slot, "ID")))
  # convert to json
  gj <- sp_to_geojson_rcpp(sp, data, use_ID, precision)
  Encoding(gj) <- "UTF-8"
  # and add class
  structure(gj,class=c("json", "geo_json"))
}