1 #include <Rcpp.h>
2 //#include "spatRaster.h"
3 #include "spatRasterMultiple.h"
4
5 //#include <memory> //std::addressof
6 #include "gdal_priv.h"
7 #include "gdalio.h"
8 #include "ogr_spatialref.h"
9
10
11 #if GDAL_VERSION_MAJOR >= 3
12 #include "proj.h"
13 #endif
14
15 // [[Rcpp::export(name = ".sameSRS")]]
sameSRS(std::string x,std::string y)16 bool sameSRS(std::string x, std::string y) {
17 std::string msg;
18 SpatSRS srs;
19 if (!srs.set(x, msg)) return false;
20 return srs.is_same(y, false);
21 }
22
23
24 // [[Rcpp::export(name = ".SRSinfo")]]
getCRSname(std::string s)25 std::vector<std::string> getCRSname(std::string s) {
26 OGRSpatialReference x;
27 OGRErr erro = x.SetFromUserInput(s.c_str());
28 if (erro != OGRERR_NONE) {
29 return {"unknown", "", "", ""};
30 }
31 std::string node;
32 if (x.IsGeographic()) {
33 node = "geogcs";
34 } else {
35 node = "projcs";
36 }
37
38 const char *value;
39 std::string name = "";
40 value = x.GetAttrValue(node.c_str());
41 if (value != NULL) {
42 name = value;
43 }
44 std::string epsg = "";
45 value = x.GetAuthorityCode(node.c_str());
46 if (value != NULL) {
47 epsg = value;
48 }
49
50 double west, south, east, north;
51 west = -10000;
52 east = -10000;
53 south = -10000;
54 north = -10000;
55
56 std::string aoi="", box="";
57 #if GDAL_VERSION_MAJOR >= 3
58 if (x.GetAreaOfUse(&west, &south, &east, &north, &value)) {
59 if (value != NULL) {
60 if (west > -1000) {
61 aoi = value;
62 box = std::to_string(west) + ", " + std::to_string(east) + ", " + std::to_string(north) + ", " + std::to_string(south);
63 }
64 }
65 }
66 #endif
67 return {name, epsg, aoi, box};
68 }
69
70 // [[Rcpp::export(name = ".getLinearUnits")]]
getLinearUnits(std::string s)71 double getLinearUnits(std::string s) {
72 std::string msg;
73 SpatSRS srs;
74 if (!srs.set(s, msg)) return NAN;
75 return srs.to_meter();
76 }
77
78
79
80 // [[Rcpp::export(name = ".geotransform")]]
geotransform(std::string fname)81 std::vector<double> geotransform(std::string fname) {
82 std::vector<double> out;
83 GDALDataset *poDataset = static_cast<GDALDataset*>(GDALOpenEx( fname.c_str(), GDAL_OF_RASTER | GDAL_OF_READONLY, NULL, NULL, NULL ));
84
85 if( poDataset == NULL ) {
86 Rcpp::Rcout << "cannot read from: " + fname << std::endl;
87 return out;
88 }
89
90 double gt[6];
91 if( poDataset->GetGeoTransform( gt ) != CE_None ) {
92 Rcpp::Rcout << "bad geotransform" << std::endl;
93 }
94 out = std::vector<double>(std::begin(gt), std::end(gt));
95 GDALClose( (GDALDatasetH) poDataset );
96
97 return out;
98 }
99
100
101 // [[Rcpp::export(name = ".gdalinfo")]]
ginfo(std::string filename,std::vector<std::string> options,std::vector<std::string> oo)102 std::string ginfo(std::string filename, std::vector<std::string> options, std::vector<std::string> oo) {
103 std::string out = gdalinfo(filename, options, oo);
104 return out;
105 }
106
107 // [[Rcpp::export(name = ".sdinfo")]]
sd_info(std::string filename)108 std::vector<std::vector<std::string>> sd_info(std::string filename) {
109 std::vector<std::vector<std::string>> sd = sdinfo(filename);
110 return sd;
111 }
112
113 // [[Rcpp::export(name = ".gdalversion")]]
gdal_version()114 std::string gdal_version() {
115 const char* what = "RELEASE_NAME";
116 const char* x = GDALVersionInfo(what);
117 std::string s = (std::string) x;
118 return s;
119 }
120
121 // [[Rcpp::export(name = ".metadata")]]
metatdata(std::string filename)122 std::vector<std::string> metatdata(std::string filename) {
123 std::vector<std::string> m = get_metadata(filename);
124 return m;
125 }
126
127 // [[Rcpp::export(name = ".sdsmetadata")]]
sdsmetatdata(std::string filename)128 std::vector<std::string> sdsmetatdata(std::string filename) {
129 std::vector<std::string> m = get_metadata_sds(filename);
130 return m;
131 }
132
133
134
135
136 // [[Rcpp::export(name = ".parsedsdsmetadata")]]
sdsmetatdataparsed(std::string filename)137 std::vector<std::vector<std::string>> sdsmetatdataparsed(std::string filename) {
138 std::vector<std::string> m = sdsmetatdata(filename);
139 std::vector<std::vector<std::string>> s = parse_metadata_sds(m);
140 return s;
141 }
142 /*
143 // [[Rcpp::export(name = ".gdaldrivers")]]
144 std::vector<std::vector<std::string>> gdal_drivers() {
145 size_t n = GetGDALDriverManager()->GetDriverCount();
146 std::vector<std::vector<std::string>> s(5);
147 for (size_t i=0; i<s.size(); i++) {
148 s[i].reserve(n);
149 }
150 GDALDriver *poDriver;
151 char **papszMetadata;
152 for (size_t i=0; i<n; i++) {
153 poDriver = GetGDALDriverManager()->GetDriver(i);
154 s[0].push_back(poDriver->GetDescription());
155 s[4].push_back(poDriver->GetMetadataItem( GDAL_DMD_LONGNAME ) );
156
157 papszMetadata = poDriver->GetMetadata();
158 bool rst = CSLFetchBoolean( papszMetadata, GDAL_DCAP_RASTER, FALSE);
159 s[1].push_back(std::to_string(rst));
160 bool create = CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE);
161 bool copy = CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE);
162 s[2].push_back(std::to_string(create + copy));
163 bool vsi = CSLFetchBoolean( papszMetadata, GDAL_DCAP_VIRTUALIO, FALSE);
164 s[3].push_back(std::to_string(vsi));
165
166 }
167 return s;
168 }
169 */
170
171
172 // [[Rcpp::export(name = ".gdaldrivers")]]
gdal_drivers()173 std::vector<std::vector<std::string>> gdal_drivers() {
174 size_t n = GetGDALDriverManager()->GetDriverCount();
175 std::vector<std::vector<std::string>> s(5, std::vector<std::string>(n));
176 GDALDriver *poDriver;
177 char **papszMetadata;
178 for (size_t i=0; i<n; i++) {
179 poDriver = GetGDALDriverManager()->GetDriver(i);
180 const char* ss = poDriver->GetDescription();
181 if (ss != NULL ) s[0][i] = ss;
182 ss = poDriver->GetMetadataItem( GDAL_DMD_LONGNAME );
183 if (ss != NULL ) s[4][i] = ss;
184
185 papszMetadata = poDriver->GetMetadata();
186 bool rst = CSLFetchBoolean( papszMetadata, GDAL_DCAP_RASTER, FALSE);
187 s[1][i] = std::to_string(rst);
188 bool create = CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE);
189 bool copy = CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE);
190 s[2][i] = std::to_string(create + copy);
191 bool vsi = CSLFetchBoolean( papszMetadata, GDAL_DCAP_VIRTUALIO, FALSE);
192 s[3][i] = std::to_string(vsi);
193 }
194 return s;
195 }
196
197
198 template <typename... Args>
warningNoCall(const char * fmt,Args &&...args)199 inline void warningNoCall(const char* fmt, Args&&... args ) {
200 Rf_warningcall(R_NilValue, tfm::format(fmt, std::forward<Args>(args)... ).c_str());
201 }
202
203 template <typename... Args>
stopNoCall(const char * fmt,Args &&...args)204 inline void NORET stopNoCall(const char* fmt, Args&&... args) {
205 throw Rcpp::exception(tfm::format(fmt, std::forward<Args>(args)... ).c_str(), false);
206 }
207
__err_warning(CPLErr eErrClass,int err_no,const char * msg)208 static void __err_warning(CPLErr eErrClass, int err_no, const char *msg) {
209 switch ( eErrClass ) {
210 case 0:
211 break;
212 case 1:
213 case 2:
214 warningNoCall("%s (GDAL %d)", msg, err_no);
215 break;
216 case 3:
217 warningNoCall("%s (GDAL error %d)", msg, err_no);
218 break;
219 case 4:
220 stopNoCall("%s (GDAL unrecoverable error %d)", msg, err_no);
221 break;
222 default:
223 warningNoCall("%s (GDAL error class %d, #%d)", msg, eErrClass, err_no);
224 break;
225 }
226 return;
227 }
228
__err_error(CPLErr eErrClass,int err_no,const char * msg)229 static void __err_error(CPLErr eErrClass, int err_no, const char *msg) {
230 switch ( eErrClass ) {
231 case 0:
232 case 1:
233 case 2:
234 break;
235 case 3:
236 warningNoCall("%s (GDAL error %d)", msg, err_no);
237 break;
238 case 4:
239 stopNoCall("%s (GDAL unrecoverable error %d)", msg, err_no);
240 break;
241 default:
242 stopNoCall("%s (GDAL unrecoverable error %d)", msg, err_no);
243 break;
244 }
245 return;
246 }
247
248
__err_fatal(CPLErr eErrClass,int err_no,const char * msg)249 static void __err_fatal(CPLErr eErrClass, int err_no, const char *msg) {
250 switch ( eErrClass ) {
251 case 0:
252 case 1:
253 case 2:
254 case 3:
255 break;
256 case 4:
257 stopNoCall("%s (GDAL unrecoverable error %d)", msg, err_no);
258 break;
259 default:
260 break;
261 }
262 return;
263 }
264
__err_none(CPLErr eErrClass,int err_no,const char * msg)265 static void __err_none(CPLErr eErrClass, int err_no, const char *msg) {
266 return;
267 }
268
269
270
271 // [[Rcpp::export(name = ".set_gdal_warnings")]]
set_gdal_warnings(int level)272 void set_gdal_warnings(int level) {
273 if (level==4) {
274 CPLSetErrorHandler((CPLErrorHandler)__err_none);
275 } else if (level==1) {
276 CPLSetErrorHandler((CPLErrorHandler)__err_warning);
277 } else if (level==2) {
278 CPLSetErrorHandler((CPLErrorHandler)__err_error);
279 } else {
280 CPLSetErrorHandler((CPLErrorHandler)__err_fatal);
281 }
282 }
283
284
285 // [[Rcpp::export(name = ".gdalinit")]]
gdal_init(std::string path)286 void gdal_init(std::string path) {
287 set_gdal_warnings(2);
288 GDALAllRegister();
289 OGRRegisterAll();
290 CPLSetConfigOption("GDAL_MAX_BAND_COUNT", "9999999");
291 //GDALregistred = true;
292 #if GDAL_VERSION_MAJOR >= 3
293 if (path != "") {
294 const char *cp = path.c_str();
295 proj_context_set_search_paths(PJ_DEFAULT_CTX, 1, &cp);
296 }
297 #endif
298 }
299
300 // [[Rcpp::export(name = ".precRank")]]
percRank(std::vector<double> x,std::vector<double> y,double minc,double maxc,int tail)301 std::vector<double> percRank(std::vector<double> x, std::vector<double> y, double minc, double maxc, int tail) {
302
303 std::vector<double> out;
304 out.reserve(y.size());
305 size_t nx = x.size();
306 for (size_t i=0; i<y.size(); i++) {
307 if (std::isnan(y[i]) ) {
308 out.push_back( NAN );
309 } else if ((y[i] < minc) | (y[i] > maxc )) {
310 out.push_back( 0 );
311 } else {
312 size_t b = 0;
313 size_t t = 0;
314 for (size_t j=0; j<x.size(); j++) {
315 if (y[i] > x[j]) {
316 b++;
317 } else if (y[i] == x[j]) {
318 t++;
319 } else {
320 // y is sorted, so we need not continue
321 break;
322 }
323 }
324 double z = (b + 0.5 * t) / nx;
325 if (tail == 1) { // both
326 if (z > 0.5) {
327 z = 2 * (1 - z);
328 } else {
329 z = 2 * z;
330 }
331 } else if (tail == 2) { // high
332 if (z < 0.5) {
333 z = 1;
334 } else {
335 z = 2 * (1 - z);
336 }
337 } else { // if (tail == 3) { // low
338 if (z > 0.5) {
339 z = 1;
340 } else {
341 z = 2 * z;
342 }
343 }
344 out.push_back(z);
345 }
346 }
347 return(out);
348 }
349
350
351 // [[Rcpp::export(name = ".setGDALCacheSizeMB")]]
setGDALCacheSizeMB(double x)352 void setGDALCacheSizeMB(double x) {
353 GDALSetCacheMax64(static_cast<int64_t>(x) * 1024 * 1024);
354 }
355
356 // [[Rcpp::export(name = ".getGDALCacheSizeMB")]]
getGDALCacheSizeMB()357 double getGDALCacheSizeMB() {
358 return static_cast<double>(GDALGetCacheMax64() / 1024 / 1024);
359 }
360
361 // convert NULL-terminated array of strings to std::vector<std::string>
charpp2vect(char ** cp)362 std::vector<std::string> charpp2vect(char **cp) {
363 std::vector<std::string> out;
364 if (cp == NULL) return out;
365 size_t i=0;
366 while (cp[i] != NULL) {
367 out.push_back(cp[i]);
368 i++;
369 }
370 return out;
371 }
372
373
374 // [[Rcpp::export(name = ".get_proj_search_paths")]]
get_proj_search_paths()375 std::vector<std::string> get_proj_search_paths() {
376 std::vector<std::string> out;
377 #if GDAL_VERSION_NUM >= 3000300
378 char **cp = OSRGetPROJSearchPaths();
379 out = charpp2vect(cp);
380 CSLDestroy(cp);
381 #else
382 out.push_back("error: GDAL >= 3.0.3 required");
383 #endif
384 return out;
385 }
386
387
388
389
390 // [[Rcpp::export(name = ".set_proj_search_paths")]]
set_proj_search_paths(std::vector<std::string> paths)391 bool set_proj_search_paths(std::vector<std::string> paths) {
392 if (!paths.size()) {
393 return false;
394 }
395 #if GDAL_VERSION_NUM >= 3000000
396 std::vector <char *> cpaths(paths.size()+1);
397 for (size_t i = 0; i < paths.size(); i++) {
398 cpaths[i] = (char *) (paths[i].c_str());
399 }
400 cpaths[cpaths.size()] = NULL;
401 OSRSetPROJSearchPaths(cpaths.data());
402 return true;
403 #else
404 return false;
405 #endif
406 }
407
408
409
410