1 #include <ogr_api.h>
2 #include <ogr_geometry.h>
3
4 #include <Rcpp.h>
5
6 #include "gdal_sf_pkg.h"
7
8 // [[Rcpp::export]]
CPL_area(Rcpp::List sfc)9 Rcpp::NumericVector CPL_area(Rcpp::List sfc) {
10 std::vector<OGRGeometry *> g = ogr_from_sfc(sfc, NULL);
11 Rcpp::NumericVector out(sfc.length());
12 for (size_t i = 0; i < g.size(); i++) {
13 if (g[i]->getDimension() == 2) {
14 OGRwkbGeometryType gt = OGR_GT_Flatten(g[i]->getGeometryType());
15 if (gt == wkbMultiSurface || gt == wkbMultiPolygon) {
16 OGRGeometryCollection *gc = (OGRGeometryCollection *) g[i];
17 out[i] = gc->get_Area();
18 } else {
19 OGRSurface *surf = (OGRSurface *) g[i];
20 out[i] = surf->get_Area();
21 }
22 } else
23 out[i] = 0.0;
24 OGRGeometryFactory::destroyGeometry(g[i]);
25 }
26 return out;
27 }
28
29 // [[Rcpp::export]]
CPL_gdal_dimension(Rcpp::List sfc,bool NA_if_empty=true)30 Rcpp::IntegerVector CPL_gdal_dimension(Rcpp::List sfc, bool NA_if_empty = true) {
31 std::vector<OGRGeometry *> g = ogr_from_sfc(sfc, NULL);
32 Rcpp::IntegerVector out(sfc.length());
33 for (size_t i = 0; i < g.size(); i++) {
34 if (NA_if_empty && g[i]->IsEmpty())
35 out[i] = NA_INTEGER;
36 else
37 out[i] = g[i]->getDimension();
38 OGRGeometryFactory f;
39 f.destroyGeometry(g[i]);
40 }
41 return out;
42 }
43
44 // [[Rcpp::export]]
CPL_length(Rcpp::List sfc)45 Rcpp::NumericVector CPL_length(Rcpp::List sfc) {
46 std::vector<OGRGeometry *> g = ogr_from_sfc(sfc, NULL);
47 Rcpp::NumericVector out(sfc.length());
48 for (size_t i = 0; i < g.size(); i++) {
49 OGRwkbGeometryType gt = OGR_GT_Flatten(g[i]->getGeometryType());
50 switch (gt) {
51 case wkbPoint:
52 case wkbMultiPoint:
53 case wkbPolygon:
54 case wkbMultiPolygon:
55 out[i] = 0.0;
56 break;
57 case wkbLineString:
58 case wkbCircularString:
59 case wkbCompoundCurve:
60 case wkbCurve: {
61 OGRCurve *a = (OGRCurve *) g[i];
62 out[i] = a->get_Length();
63 }
64 break;
65 default: {
66 OGRGeometryCollection *a = (OGRGeometryCollection *) g[i];
67 out[i] = a->get_Length();
68 }
69 }
70 OGRGeometryFactory f;
71 f.destroyGeometry(g[i]);
72 }
73 return out;
74 }
75
76 // [[Rcpp::export]]
CPL_gdal_segmentize(Rcpp::List sfc,double dfMaxLength=0.0)77 Rcpp::List CPL_gdal_segmentize(Rcpp::List sfc, double dfMaxLength = 0.0) {
78
79 if (dfMaxLength <= 0.0)
80 Rcpp::stop("argument dfMaxLength should be positive\n");
81
82 std::vector<OGRGeometry *> g = ogr_from_sfc(sfc, NULL);
83 for (size_t i = 0; i < g.size(); i++)
84 g[i]->segmentize(dfMaxLength);
85 Rcpp::List ret = sfc_from_ogr(g, true);
86 ret.attr("crs") = sfc.attr("crs");
87 return ret;
88 }
89
90 // [[Rcpp::export]]
CPL_gdal_linestring_sample(Rcpp::List sfc,Rcpp::List distLst)91 Rcpp::List CPL_gdal_linestring_sample(Rcpp::List sfc, Rcpp::List distLst) {
92 if (sfc.size() != distLst.size())
93 Rcpp::stop("sfc and dist should have equal length"); // #nocov
94 std::vector<OGRGeometry *> g = ogr_from_sfc(sfc, NULL);
95 std::vector<OGRGeometry *> out(g.size());
96 for (size_t i = 0; i < g.size(); i++) {
97 if (wkbFlatten(g[i]->getGeometryType()) != wkbLineString)
98 Rcpp::stop("CPL_gdal_linestring_sample only available for LINESTRING"); // #nocov
99 OGRGeometryCollection *gc = new OGRGeometryCollection;
100 Rcpp::NumericVector dists = distLst[i];
101 for (int j = 0; j < dists.size(); j++) {
102 OGRPoint *poPoint = new OGRPoint;
103 ((OGRLineString *) g[i])->Value(dists[j], poPoint);
104 gc->addGeometryDirectly(poPoint);
105 }
106 out[i] = OGRGeometryFactory::forceToMultiPoint(gc);
107 }
108 Rcpp::List ret = sfc_from_ogr(g, true); // releases g
109 ret = sfc_from_ogr(out, true); // releases out
110 ret.attr("crs") = sfc.attr("crs");
111 return ret;
112 }
113