1 /* Copyright (c) 2019-20 Roger Bivand */
2
3
4 #include <R.h>
5 #include <Rdefines.h>
6 #include "rgdal.h"
7 //#ifdef ACCEPT_USE_OF_DEPRECATED_PROJ_API_H // kludge for 6 only
8 #ifdef PROJ_H_API // PROJ_H_API
9 #include <proj.h>
10
11 #ifdef __cplusplus
12 extern "C" {
13 #endif
14
silent_logger(void *,int,const char *)15 static void silent_logger(void *, int, const char *) {}
16
17
RGDAL_checkCRSArgs(SEXP args)18 SEXP RGDAL_checkCRSArgs(SEXP args) {
19 // Rprintf("Not available for PROJ version >= 6");
20 return(R_NilValue);
21 }
22
PROJ4NADsInstalled(void)23 SEXP PROJ4NADsInstalled(void) {
24 // Rprintf("Not available for PROJ version >= 6");
25 return(R_NilValue);
26 }
27
PROJ4_proj_def_dat_Installed(void)28 SEXP PROJ4_proj_def_dat_Installed(void) {
29 // Rprintf("Not available for PROJ version >= 6");
30 return(R_NilValue);
31 }
32
transform(SEXP fromargs,SEXP toargs,SEXP npts,SEXP x,SEXP y,SEXP z)33 SEXP transform(SEXP fromargs, SEXP toargs, SEXP npts, SEXP x, SEXP y, SEXP z) {
34 // Rprintf("Not available for PROJ version >= 6");
35 return(R_NilValue);
36 }
37
38
RGDAL_project(SEXP n,SEXP xlon,SEXP ylat,SEXP projarg,SEXP ob_tran)39 SEXP RGDAL_project(SEXP n, SEXP xlon, SEXP ylat, SEXP projarg, SEXP ob_tran) {
40 // Rprintf("Not available for PROJ version >= 6");
41 return(R_NilValue);
42 }
43
44
project_inv(SEXP n,SEXP x,SEXP y,SEXP projarg,SEXP ob_tran)45 SEXP project_inv(SEXP n, SEXP x, SEXP y, SEXP projarg, SEXP ob_tran) {
46 // Rprintf("Not available for PROJ version >= 6");
47 return(R_NilValue);
48 }
49
50
51
PROJ4VersionInfo(void)52 SEXP PROJ4VersionInfo(void) {
53 SEXP ans;
54 PJ_INFO pji= proj_info();
55
56 PROTECT(ans=NEW_LIST(2));
57 SET_VECTOR_ELT(ans, 0, NEW_CHARACTER(1));
58 SET_VECTOR_ELT(ans, 1, NEW_INTEGER(1));
59 SET_STRING_ELT(VECTOR_ELT(ans, 0), 0,
60 COPY_TO_USER_STRING(pji.release));
61 INTEGER_POINTER(VECTOR_ELT(ans, 1))[0] = (pji.major*100)+(pji.minor*10)+pji.patch;
62
63 UNPROTECT(1);
64
65 return(ans);
66 }
67
RGDAL_projInfo(SEXP type)68 SEXP RGDAL_projInfo(SEXP type) {
69 SEXP ans=NULL;
70 SEXP ansnames;
71 int n=0, pc=0;
72
73
74 if (INTEGER_POINTER(type)[0] == 0) {
75 PROTECT(ans = NEW_LIST(2)); pc++;
76 PROTECT(ansnames = NEW_CHARACTER(2)); pc++;
77 SET_STRING_ELT(ansnames, 0, COPY_TO_USER_STRING("name"));
78 SET_STRING_ELT(ansnames, 1, COPY_TO_USER_STRING("description"));
79 setAttrib(ans, R_NamesSymbol, ansnames);
80
81 const struct PJ_LIST *lp;
82 for (lp = proj_list_operations() ; lp->id ; ++lp) {
83 if( strcmp(lp->id,"latlong") == 0
84 || strcmp(lp->id,"longlat") == 0
85 || strcmp(lp->id,"geocent") == 0 )
86 continue;
87 n++;
88 }
89 SET_VECTOR_ELT(ans, 0, NEW_CHARACTER(n));
90 SET_VECTOR_ELT(ans, 1, NEW_CHARACTER(n));
91 n=0;
92 for (lp = proj_list_operations() ; lp->id ; ++lp) {
93 if( strcmp(lp->id,"latlong") == 0
94 || strcmp(lp->id,"longlat") == 0
95 || strcmp(lp->id,"geocent") == 0 )
96 continue;
97 SET_STRING_ELT(VECTOR_ELT(ans, 0), n,
98 COPY_TO_USER_STRING(lp->id));
99
100 SET_STRING_ELT(VECTOR_ELT(ans, 1), n,
101 COPY_TO_USER_STRING(*lp->descr));
102 n++;
103 }
104 } else if (INTEGER_POINTER(type)[0] == 1) {
105 PROTECT(ans = NEW_LIST(4)); pc++;
106 PROTECT(ansnames = NEW_CHARACTER(4)); pc++;
107 SET_STRING_ELT(ansnames, 0, COPY_TO_USER_STRING("name"));
108 SET_STRING_ELT(ansnames, 1, COPY_TO_USER_STRING("major"));
109 SET_STRING_ELT(ansnames, 2, COPY_TO_USER_STRING("ell"));
110 SET_STRING_ELT(ansnames, 3, COPY_TO_USER_STRING("description"));
111 setAttrib(ans, R_NamesSymbol, ansnames);
112
113 const struct PJ_ELLPS *le;
114 for (le = proj_list_ellps(); le->id ; ++le) n++;
115 SET_VECTOR_ELT(ans, 0, NEW_CHARACTER(n));
116 SET_VECTOR_ELT(ans, 1, NEW_CHARACTER(n));
117 SET_VECTOR_ELT(ans, 2, NEW_CHARACTER(n));
118 SET_VECTOR_ELT(ans, 3, NEW_CHARACTER(n));
119 n=0;
120 for (le = proj_list_ellps(); le->id ; ++le) {
121 SET_STRING_ELT(VECTOR_ELT(ans, 0), n,
122 COPY_TO_USER_STRING(le->id));
123 SET_STRING_ELT(VECTOR_ELT(ans, 1), n,
124 COPY_TO_USER_STRING(le->major));
125 SET_STRING_ELT(VECTOR_ELT(ans, 2), n,
126 COPY_TO_USER_STRING(le->ell));
127 SET_STRING_ELT(VECTOR_ELT(ans, 3), n,
128 COPY_TO_USER_STRING(le->name));
129 n++;
130 }
131 } else if (INTEGER_POINTER(type)[0] == 2) {
132 return(R_NilValue);
133 } else if (INTEGER_POINTER(type)[0] == 3) {
134 PROTECT(ans = NEW_LIST(3)); pc++;
135 PROTECT(ansnames = NEW_CHARACTER(3)); pc++;
136 SET_STRING_ELT(ansnames, 0, COPY_TO_USER_STRING("id"));
137 SET_STRING_ELT(ansnames, 1, COPY_TO_USER_STRING("to_meter"));
138 SET_STRING_ELT(ansnames, 2, COPY_TO_USER_STRING("name"));
139 setAttrib(ans, R_NamesSymbol, ansnames);
140 #if ((PROJ_VERSION_MAJOR == 7 && PROJ_VERSION_MINOR >= 1) || PROJ_VERSION_MAJOR > 7)
141 PROJ_UNIT_INFO** units;
142 units = proj_get_units_from_database(nullptr, nullptr, "linear", false, nullptr);
143 n = 0;
144 //Rprintf("n: %d\n", n);
145 for (int i = 0; units && units[i]; i++) {
146 if (units[i]->proj_short_name) n++;
147 //Rprintf("%s\n", units[i]->proj_short_name);
148 }
149 //Rprintf("n: %d\n", n);
150 SET_VECTOR_ELT(ans, 0, NEW_CHARACTER(n));
151 SET_VECTOR_ELT(ans, 1, NEW_NUMERIC(n));
152 SET_VECTOR_ELT(ans, 2, NEW_CHARACTER(n));
153 int n1=0;
154 for (int i = 0; units && units[i]; i++) {
155 if (units[i]->proj_short_name) {
156 SET_STRING_ELT(VECTOR_ELT(ans, 0), n1,
157 COPY_TO_USER_STRING(units[i]->proj_short_name));
158 NUMERIC_POINTER(VECTOR_ELT(ans, 1))[n1] = units[i]->conv_factor;
159 SET_STRING_ELT(VECTOR_ELT(ans, 2), n1,
160 COPY_TO_USER_STRING(units[i]->name));
161 n1++;
162 //Rprintf("%s %f %s\n", units[i]->proj_short_name, units[i]->conv_factor, units[i]->name);
163 }
164 if (n1 >= n) break; //EPJ correction
165 }
166 proj_unit_list_destroy(units);
167 // >= 710 proj_get_units_from_database()
168 #else
169 const struct PJ_UNITS *lu;
170 for (lu = proj_list_units(); lu->id ; ++lu) n++;
171 SET_VECTOR_ELT(ans, 0, NEW_CHARACTER(n));
172 SET_VECTOR_ELT(ans, 1, NEW_CHARACTER(n));
173 SET_VECTOR_ELT(ans, 2, NEW_CHARACTER(n));
174 n=0;
175 for (lu = proj_list_units(); lu->id ; ++lu) {
176 SET_STRING_ELT(VECTOR_ELT(ans, 0), n,
177 COPY_TO_USER_STRING(lu->id));
178 SET_STRING_ELT(VECTOR_ELT(ans, 1), n,
179 COPY_TO_USER_STRING(lu->to_meter));
180 SET_STRING_ELT(VECTOR_ELT(ans, 2), n,
181 COPY_TO_USER_STRING(lu->name));
182 n++;
183 }
184 #endif
185 } else error("no such type");
186
187 UNPROTECT(pc);
188 return(ans);
189 }
190
get_proj_search_path(void)191 SEXP get_proj_search_path(void) {
192
193 SEXP res;
194 PROTECT(res = NEW_CHARACTER(1));
195 SET_STRING_ELT(res, 0, COPY_TO_USER_STRING(proj_info().searchpath));
196 UNPROTECT(1);
197 return(res);
198
199 }
200
get_proj_user_writable_dir()201 SEXP get_proj_user_writable_dir() {
202
203 #if ((PROJ_VERSION_MAJOR == 7 && PROJ_VERSION_MINOR >= 1) || PROJ_VERSION_MAJOR > 7)
204 SEXP res;
205 PROTECT(res = NEW_CHARACTER(1));
206 SET_STRING_ELT(res, 0, COPY_TO_USER_STRING(proj_context_get_user_writable_directory(PJ_DEFAULT_CTX, false)));
207 UNPROTECT(1);
208 return(res);
209 #else
210 return(R_NilValue);
211 #endif
212
213 }
214
215
set_proj_paths(SEXP paths)216 SEXP set_proj_paths(SEXP paths) {
217 //PJ_CONTEXT *ctx = proj_context_create();
218 SEXP res;
219 int i, n = length(paths);
220 char **paths_loc = (char **) R_alloc((size_t) n, sizeof(char*));
221 for (i=0; i<n; i++) paths_loc[i] = (char *) CHAR(STRING_ELT(paths, i));
222 // const char *paths = CHAR(STRING_ELT(path, 0));
223 proj_context_set_search_paths(PJ_DEFAULT_CTX, n, (const char* const*) paths_loc);
224 if (int this_errno = proj_context_errno(PJ_DEFAULT_CTX) != 0) {
225 const char *errstr = proj_errno_string(this_errno);
226 //proj_context_destroy(ctx);
227 error("setting search path failed: %s", errstr);
228 }
229 PROTECT(res = NEW_CHARACTER(1));
230 SET_STRING_ELT(res, 0, COPY_TO_USER_STRING(proj_info().searchpath));
231 UNPROTECT(1);
232 //proj_context_destroy(ctx);
233 return(res);
234 }
235
get_source_crs(SEXP source)236 SEXP get_source_crs(SEXP source) {
237
238 PJ_CONTEXT *ctx = proj_context_create();
239 PJ *source_crs, *target_crs;
240 SEXP res;
241
242 source_crs = proj_create(ctx, CHAR(STRING_ELT(source, 0)));
243
244 if (source_crs == NULL) {
245 proj_context_destroy(ctx);
246 error("source crs not created");
247 }
248
249 target_crs = proj_get_source_crs(ctx, source_crs);
250
251 if (target_crs == NULL) {
252 proj_context_destroy(ctx);
253 error("target crs not created");
254 }
255
256 PROTECT(res = NEW_CHARACTER(1));
257 SET_STRING_ELT(res, 0,
258 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 3
259 COPY_TO_USER_STRING(proj_as_wkt(ctx, target_crs, PJ_WKT2_2018, NULL))
260 #else
261 COPY_TO_USER_STRING(proj_as_wkt(ctx, target_crs, PJ_WKT2_2019, NULL))
262 #endif
263 );
264 UNPROTECT(1);
265 proj_destroy(target_crs);
266 proj_destroy(source_crs);
267 proj_context_destroy(ctx);
268
269 return(res);
270
271
272 }
273
proj_vis_order(SEXP wkt2)274 SEXP proj_vis_order(SEXP wkt2) {
275
276 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 3
277 warning("no CRS normalization available before PROJ 6.3");
278 return(wkt2);
279 #else
280
281 PJ_CONTEXT *ctx = proj_context_create();
282 PJ *source_crs, *target_crs;
283 SEXP res;
284
285 source_crs = proj_create(ctx, CHAR(STRING_ELT(wkt2, 0)));
286
287 if (source_crs == NULL) {
288 proj_context_destroy(ctx);
289 error("proj_vis_order: source crs not created");
290 }
291 target_crs = proj_normalize_for_visualization(ctx, source_crs);
292 if (target_crs == NULL) {
293 proj_context_destroy(ctx);
294 error("proj_vis_order: target crs not created");
295 }
296
297 PROTECT(res = NEW_CHARACTER(1));
298 SET_STRING_ELT(res, 0,
299 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 3
300 COPY_TO_USER_STRING(proj_as_wkt(ctx, target_crs, PJ_WKT2_2018, NULL))
301 #else
302 COPY_TO_USER_STRING(proj_as_wkt(ctx, target_crs, PJ_WKT2_2019, NULL))
303 #endif
304 );
305 UNPROTECT(1);
306 proj_destroy(target_crs);
307 proj_destroy(source_crs);
308 proj_context_destroy(ctx);
309
310 return(res);
311 #endif
312 }
313
314 // unname(sapply(o[[2]], function(x) gsub(" ", " +", paste0("+", x))))
315
list_coordinate_ops(SEXP source,SEXP target,SEXP area_of_interest,SEXP strict_containment,SEXP vis_order)316 SEXP list_coordinate_ops(SEXP source, SEXP target, SEXP area_of_interest, SEXP strict_containment, SEXP vis_order) {
317
318 PJ_CONTEXT *ctx = proj_context_create();
319 PJ_OPERATION_FACTORY_CONTEXT* operation_factory_context = NULL;
320 PJ_OBJ_LIST *pj_operations = NULL;
321 PJ *source_crs, *target_crs;
322 // PJ* pj_transform = NULL;
323 SEXP ans, input;
324 int num_operations, i, j, is_instantiable, is_ballpark, grid_count;
325 int pc=0;
326 double accuracy;
327 int grid_OK, out_direct_download, out_open_license, out_available;
328 const char *out_short_name, *out_full_name, *out_package_name, *out_url;
329 PJ_PROJ_INFO pjinfo;
330
331 // valgrind 210120
332 source_crs = proj_create(ctx, CHAR(STRING_ELT(source, 0)));
333
334 if (source_crs == NULL) {
335 proj_context_destroy(ctx);
336 error("source crs not created");
337 }
338
339 // valgrind 210120
340 target_crs = proj_create(ctx, CHAR(STRING_ELT(target, 0)));
341
342 if (target_crs == NULL) {
343 proj_destroy(source_crs);
344 proj_context_destroy(ctx);
345 error("target crs not created");
346 }
347
348 operation_factory_context = proj_create_operation_factory_context(ctx,
349 NULL);
350 if (operation_factory_context == NULL) {
351 proj_destroy(source_crs); proj_destroy(target_crs);
352 proj_context_destroy(ctx);
353 error("operation factory context not created");
354 }
355
356 if (!ISNA(NUMERIC_POINTER(area_of_interest)[0])) {
357 proj_operation_factory_context_set_area_of_interest(
358 ctx, operation_factory_context,
359 NUMERIC_POINTER(area_of_interest)[0],
360 NUMERIC_POINTER(area_of_interest)[1],
361 NUMERIC_POINTER(area_of_interest)[2],
362 NUMERIC_POINTER(area_of_interest)[3]);
363 }
364
365 if (LOGICAL_POINTER(strict_containment)[0])
366 proj_operation_factory_context_set_spatial_criterion(
367 ctx, operation_factory_context,
368 PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT);
369 else proj_operation_factory_context_set_spatial_criterion(
370 ctx, operation_factory_context,
371 PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
372
373 proj_operation_factory_context_set_grid_availability_use(
374 ctx, operation_factory_context,
375 PROJ_GRID_AVAILABILITY_USED_FOR_SORTING);
376
377 // valgrind 210120
378 pj_operations = proj_create_operations(ctx,
379 source_crs, target_crs, operation_factory_context);
380
381 if (pj_operations == NULL) {
382 proj_operation_factory_context_destroy(operation_factory_context);
383 proj_destroy(source_crs); proj_destroy(target_crs);
384 proj_context_destroy(ctx);
385 error("operations list not created");
386 }
387
388 num_operations = proj_list_get_count(pj_operations);
389
390 if (num_operations < 1L) {
391 proj_list_destroy(pj_operations);
392 proj_operation_factory_context_destroy(operation_factory_context);
393 proj_destroy(source_crs); proj_destroy(target_crs);
394 proj_context_destroy(ctx);
395 return(R_NilValue);
396 }
397
398 PROTECT(ans=NEW_LIST(7)); pc++;
399 SET_VECTOR_ELT(ans, 0, NEW_CHARACTER(num_operations));
400 SET_VECTOR_ELT(ans, 1, NEW_CHARACTER(num_operations));
401 SET_VECTOR_ELT(ans, 2, NEW_NUMERIC(num_operations));
402 SET_VECTOR_ELT(ans, 3, NEW_LOGICAL(num_operations));
403 SET_VECTOR_ELT(ans, 4, NEW_LOGICAL(num_operations));
404 SET_VECTOR_ELT(ans, 5, NEW_INTEGER(num_operations));
405 SET_VECTOR_ELT(ans, 6, NEW_LIST(num_operations));
406
407 PROTECT(input=NEW_LIST(5)); pc++;
408 SET_VECTOR_ELT(input, 0, source);
409 SET_VECTOR_ELT(input, 1, target);
410 SET_VECTOR_ELT(input, 2, area_of_interest);
411 SET_VECTOR_ELT(input, 3, strict_containment);
412 SET_VECTOR_ELT(input, 4, vis_order);
413 setAttrib(ans, install("input"), input);
414
415
416 for (i=0; i<num_operations; i++) {
417 PJ * pj_transform = proj_list_get(ctx, pj_operations, i);
418 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 1
419 warning("no Coordinate Operation normalization available before PROJ 6.1");
420 #else
421 if (LOGICAL_POINTER(vis_order)[0]) {
422 PJ* pj_trans_orig = pj_transform;
423 pj_transform = proj_normalize_for_visualization(ctx, pj_trans_orig);
424 proj_destroy(pj_trans_orig);
425 }
426 #endif
427 is_instantiable = proj_coordoperation_is_instantiable(ctx,
428 pj_transform);
429 is_ballpark = proj_coordoperation_has_ballpark_transformation(ctx,
430 pj_transform);
431 accuracy = proj_coordoperation_get_accuracy(ctx,
432 pj_transform);
433 grid_count = proj_coordoperation_get_grid_used_count(ctx,
434 pj_transform);
435 pjinfo = proj_pj_info(pj_transform);
436
437 SET_STRING_ELT(VECTOR_ELT(ans, 0), i,
438 COPY_TO_USER_STRING(pjinfo.description));
439 SET_STRING_ELT(VECTOR_ELT(ans, 1), i,
440 COPY_TO_USER_STRING(pjinfo.definition));
441
442 NUMERIC_POINTER(VECTOR_ELT(ans, 2))[i] = accuracy;
443 LOGICAL_POINTER(VECTOR_ELT(ans, 3))[i] = is_instantiable;
444 LOGICAL_POINTER(VECTOR_ELT(ans, 4))[i] = is_ballpark;
445 INTEGER_POINTER(VECTOR_ELT(ans, 5))[i] = grid_count;
446
447 if (grid_count > 0L) {
448 SET_VECTOR_ELT(VECTOR_ELT(ans, 6), i, NEW_LIST(grid_count));
449 for (j=0; j<grid_count; j++) {
450 grid_OK = proj_coordoperation_get_grid_used(ctx, pj_transform,
451 j, &out_short_name, &out_full_name, &out_package_name,
452 &out_url, &out_direct_download, &out_open_license,
453 &out_available);
454 if (grid_OK) {
455 SET_VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(ans, 6), i), j,
456 NEW_LIST(7));
457 SET_VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT( ans, 6),
458 i), j), 0, NEW_CHARACTER(1));
459 SET_STRING_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(
460 VECTOR_ELT(ans, 6), i), j), 0), 0,
461 COPY_TO_USER_STRING(out_short_name));
462 SET_VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT( ans, 6),
463 i), j), 1, NEW_CHARACTER(1));
464 SET_STRING_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(
465 VECTOR_ELT(ans, 6), i), j), 1), 0,
466 COPY_TO_USER_STRING(out_full_name));
467 SET_VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT( ans, 6),
468 i), j), 2, NEW_CHARACTER(1));
469 SET_STRING_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(
470 VECTOR_ELT(ans, 6), i), j), 2), 0,
471 COPY_TO_USER_STRING(out_package_name));
472 SET_VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT( ans, 6),
473 i), j), 3, NEW_CHARACTER(1));
474 SET_STRING_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(
475 VECTOR_ELT(ans, 6), i), j), 3), 0,
476 COPY_TO_USER_STRING(out_url));
477 SET_VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT( ans, 6),
478 i), j), 4, NEW_LOGICAL(1));
479 LOGICAL_POINTER(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(
480 VECTOR_ELT(ans, 6), i), j), 4))[0] =
481 out_direct_download;
482 SET_VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT( ans, 6),
483 i), j), 5, NEW_LOGICAL(1));
484 LOGICAL_POINTER(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(
485 VECTOR_ELT(ans, 6), i), j), 5))[0] =
486 out_open_license;
487 SET_VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT( ans, 6),
488 i), j), 6, NEW_LOGICAL(1));
489 LOGICAL_POINTER(VECTOR_ELT(VECTOR_ELT(VECTOR_ELT(
490 VECTOR_ELT(ans, 6), i), j), 6))[0] =
491 out_available;
492
493 }
494 }
495 }
496
497 proj_destroy(pj_transform);
498
499 }
500
501 // proj_destroy(pj_transform);
502 proj_list_destroy(pj_operations);
503 proj_operation_factory_context_destroy(operation_factory_context);
504 proj_destroy(source_crs); proj_destroy(target_crs);
505 proj_context_destroy(ctx);
506
507 UNPROTECT(pc);
508 return(ans);
509
510 }
511
CRS_compare(SEXP fromargs,SEXP toargs,SEXP type1,SEXP type2)512 SEXP CRS_compare(SEXP fromargs, SEXP toargs, SEXP type1, SEXP type2) {
513
514 //PJ_CONTEXT *ctx = proj_context_create();
515 PJ *source_crs, *target_crs;
516 SEXP res;
517 int ires_strict, ires_equiv, ires_equiv_ao;
518 //Rprintf("source crs input: %s\n", CHAR(STRING_ELT(fromargs, 0)));
519 if ((source_crs = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(fromargs, 0)))) == NULL) {
520 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
521 //proj_context_destroy(ctx);
522 error("source crs creation failed: %s", errstr);
523 }
524
525 //Rprintf("source crs: %s\n", proj_as_proj_string(PJ_DEFAULT_CTX, source_crs, PJ_PROJ_5, NULL)); not filled for WKT
526 //Rprintf("target crs input: %s\n", CHAR(STRING_ELT(toargs, 0)));
527 if ((target_crs = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(toargs, 0)))) == NULL) {
528 proj_destroy(source_crs);
529 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
530 //proj_context_destroy(ctx);
531 error("target crs creation failed: %s", errstr);
532 }
533 //Rprintf("target crs: %s\n", proj_as_proj_string(PJ_DEFAULT_CTX, target_crs, PJ_PROJ_5, NULL)); not filled for WKT
534 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 3
535 ires_strict = proj_is_equivalent_to(source_crs, target_crs,
536 PJ_COMP_STRICT);
537 ires_equiv = proj_is_equivalent_to(source_crs, target_crs,
538 PJ_COMP_EQUIVALENT);
539 ires_equiv_ao = proj_is_equivalent_to(source_crs, target_crs,
540 PJ_COMP_EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS);
541 #else
542 ires_strict = proj_is_equivalent_to_with_ctx(PJ_DEFAULT_CTX, source_crs, target_crs,
543 PJ_COMP_STRICT);
544 ires_equiv = proj_is_equivalent_to_with_ctx(PJ_DEFAULT_CTX, source_crs, target_crs,
545 PJ_COMP_EQUIVALENT);
546 ires_equiv_ao = proj_is_equivalent_to_with_ctx(PJ_DEFAULT_CTX, source_crs, target_crs,
547 PJ_COMP_EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS);
548 #endif
549 //Rprintf("ires: %d, iresao\n", ires);
550 PROTECT(res = NEW_INTEGER(3));
551 INTEGER_POINTER(res)[0] = ires_strict;
552 INTEGER_POINTER(res)[1] = ires_equiv;
553 INTEGER_POINTER(res)[2] = ires_equiv_ao;
554 proj_destroy(target_crs);
555 proj_destroy(source_crs);
556 //proj_context_destroy(ctx);
557 UNPROTECT(1);
558 return(res);
559 }
560
transform_ng(SEXP fromargs,SEXP toargs,SEXP coordOp,SEXP npts,SEXP x,SEXP y,SEXP z,SEXP aoi)561 SEXP transform_ng(SEXP fromargs, SEXP toargs, SEXP coordOp, SEXP npts, SEXP x, SEXP y, SEXP z, SEXP aoi) {
562
563 //PJ_CONTEXT *ctx = proj_context_create();
564 PJ *source_crs = NULL, *target_crs = NULL;
565 PJ* pj_transform = NULL;
566 PJ_AREA *area_of_interest = 0;
567
568 int i, n, nwarn=0, have_z, pc=0, have_CO, vis_order, use_aoi=1;
569 double *xx, *yy, *zz=NULL;
570 SEXP enforce_xy = getAttrib(npts, install("enforce_xy"));
571 SEXP res;
572
573 if (coordOp == R_NilValue) have_CO = 0;
574 else have_CO = 1;
575
576 if (z == R_NilValue) have_z = 0;
577 else have_z = 1;
578
579
580 if (enforce_xy == R_NilValue) vis_order = 0;
581 else if (LOGICAL_POINTER(enforce_xy)[0] == 1) vis_order = 1;
582 else if (LOGICAL_POINTER(enforce_xy)[0] == 0) vis_order = 0;
583 else vis_order = 0;
584
585 if (aoi == R_NilValue) use_aoi = 0;
586
587 if (!have_CO && use_aoi) {
588 area_of_interest = proj_area_create();
589 proj_area_set_bbox(area_of_interest, NUMERIC_POINTER(aoi)[0],
590 NUMERIC_POINTER(aoi)[1], NUMERIC_POINTER(aoi)[2],
591 NUMERIC_POINTER(aoi)[3]);
592 }
593 // proj_area_destroy(area_of_interest);
594
595
596
597 //Rprintf("have_CO: %d, have_z: %d: %d\n", have_CO, have_z);
598
599 if (have_CO) {
600 //Rprintf("coordinate operation input: %s\n", CHAR(STRING_ELT(coordOp, 0)));
601 if ((pj_transform = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(coordOp, 0)))) == 0) {
602 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
603 //proj_context_destroy(ctx);
604 error("coordinate operation creation failed: %s", errstr);
605 }
606 } else {
607 //Rprintf("source crs input: %s\n", CHAR(STRING_ELT(fromargs, 0)));
608 if ((source_crs = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(fromargs, 0)))) == NULL) {
609 proj_area_destroy(area_of_interest);
610 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
611 //proj_context_destroy(ctx);
612 error("source crs creation failed: %s", errstr);
613 }
614
615 //Rprintf("source crs: %s\n", proj_as_proj_string(PJ_DEFAULT_CTX, source_crs, PJ_PROJ_5, NULL));
616 //Rprintf("target crs input: %s\n", CHAR(STRING_ELT(toargs, 0)));
617 if ((target_crs = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(toargs, 0)))) == NULL) {
618 proj_area_destroy(area_of_interest);
619 proj_destroy(source_crs);
620 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
621 //proj_context_destroy(ctx);
622 error("target crs creation failed: %s", errstr);
623 }
624 //Rprintf("target crs: %s\n", proj_as_proj_string(PJ_DEFAULT_CTX, target_crs, PJ_PROJ_5, NULL));
625 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 2
626 if ((pj_transform = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
627 CHAR(STRING_ELT(fromargs, 0)),
628 CHAR(STRING_ELT(toargs, 0)), area_of_interest)) == 0) {
629 // FIXME >= 6.2.0
630 proj_area_destroy(area_of_interest);
631 proj_destroy(target_crs);
632 proj_destroy(source_crs);
633 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
634 //proj_context_destroy(ctx);
635 error("coordinate operation creation from WKT failed: %s", errstr);
636 }
637
638 #else
639 if ((pj_transform = proj_create_crs_to_crs_from_pj(PJ_DEFAULT_CTX,
640 source_crs, target_crs, area_of_interest, NULL)) == 0) {
641 // FIXME >= 6.2.0
642 proj_area_destroy(area_of_interest);
643 proj_destroy(target_crs);
644 proj_destroy(source_crs);
645 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
646 //proj_context_destroy(ctx);
647 error("coordinate operation creation from WKT failed: %s", errstr);
648 }
649 #endif
650 //Rprintf("%s\n", proj_pj_info(pj_transform).definition);
651 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 1
652 warning("no Coordinate Operation normalization available before PROJ 6.1");
653 #else
654 if (vis_order == 1) {
655 PJ* pj_trans_orig = pj_transform;
656 pj_transform = proj_normalize_for_visualization(PJ_DEFAULT_CTX, pj_trans_orig);
657 proj_destroy(pj_trans_orig);
658 }
659 #endif
660 }
661 //Rprintf("%s\n", proj_pj_info(pj_transform).definition);
662 // if (proj_normalize_for_visualization(ctx, pj_transform) != NULL) // EJP
663 // pj_transform = proj_normalize_for_visualization(ctx, pj_transform); // EJP
664 //Rprintf("%s\n", proj_pj_info(pj_transform).definition);
665
666 n = INTEGER_POINTER(npts)[0];
667 xx = (double *) R_alloc((size_t) n, sizeof(double));
668 yy = (double *) R_alloc((size_t) n, sizeof(double));
669 if (have_z) zz = (double *) R_alloc((size_t) n, sizeof(double));
670
671 for (i=0; i < n; i++) {
672 xx[i] = NUMERIC_POINTER(x)[i];
673 yy[i] = NUMERIC_POINTER(y)[i];
674 if (have_z) zz[i] = NUMERIC_POINTER(z)[i];
675 }
676 /* if (ob_tran == 1) {
677 for (i=0; i < n; i++) {
678 xx[i] = proj_torad(xx[i]);
679 yy[i] = proj_torad(yy[i]);
680 }
681 }*/
682 size_t stride = sizeof(double), n1;
683 // FIXME handle out-of-domain output https://github.com/r-spatial/sf/issues/1223
684 if (have_z) {
685 if((n1 = proj_trans_generic(pj_transform, PJ_FWD, xx, stride,
686 (size_t) n, yy, stride, (size_t) n, zz, stride, (size_t) n,
687 NULL, stride, 0)) != (size_t) n) {
688 proj_destroy(pj_transform);
689 if (!have_CO) {
690 proj_area_destroy(area_of_interest);
691 proj_destroy(target_crs);
692 proj_destroy(source_crs);
693 }
694 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
695 //proj_context_destroy(ctx);
696 error("error in proj_transform_generic: %n of %n coordinates succeeded\n %s", n1, n, errstr);
697 }
698 } else {
699 if((n1 = proj_trans_generic(pj_transform, PJ_FWD, xx, stride,
700 (size_t) n, yy, stride, (size_t) n, NULL, stride, 0,
701 NULL, stride, 0)) != (size_t) n) {
702 proj_destroy(pj_transform);
703 if (!have_CO) {
704 proj_area_destroy(area_of_interest);
705 proj_destroy(target_crs);
706 proj_destroy(source_crs);
707 }
708 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
709 //proj_context_destroy(ctx);
710 error("error in proj_transform_generic: %n of %n coordinates succeeded\n %s", n1, n, errstr);
711 }
712 }
713 //Rprintf("%s\n", proj_pj_info(pj_transform).definition);
714 if (have_z) {PROTECT(res = NEW_LIST(6)); pc++;}
715 else {PROTECT(res = NEW_LIST(5)); pc++;}
716 SET_VECTOR_ELT(res, 0, NEW_NUMERIC(n));
717 SET_VECTOR_ELT(res, 1, NEW_NUMERIC(n));
718 if (have_z) {
719 SET_VECTOR_ELT(res, 2, NEW_NUMERIC(n));
720 SET_VECTOR_ELT(res, 5, NEW_CHARACTER(1));
721 } else {
722 SET_VECTOR_ELT(res, 4, NEW_CHARACTER(1));
723 }
724
725 if (have_z) {
726 if (have_CO) {
727 SET_VECTOR_ELT(res, 3, fromargs);
728 SET_VECTOR_ELT(res, 4, toargs);
729 } else {
730 SET_VECTOR_ELT(res, 3, NEW_CHARACTER(1));
731 SET_VECTOR_ELT(res, 4, NEW_CHARACTER(1));
732 SET_STRING_ELT(VECTOR_ELT(res, 3), 0,
733 COPY_TO_USER_STRING(proj_pj_info(source_crs).definition));
734 SET_STRING_ELT(VECTOR_ELT(res, 4), 0,
735 COPY_TO_USER_STRING(proj_pj_info(target_crs).definition));
736 }
737 SET_STRING_ELT(VECTOR_ELT(res, 5), 0,
738 COPY_TO_USER_STRING(proj_pj_info(pj_transform).definition));
739 } else {
740 if (have_CO) {
741 SET_VECTOR_ELT(res, 2, fromargs);
742 SET_VECTOR_ELT(res, 3, toargs);
743 } else {
744 SET_VECTOR_ELT(res, 2, NEW_CHARACTER(1));
745 SET_VECTOR_ELT(res, 3, NEW_CHARACTER(1));
746 SET_STRING_ELT(VECTOR_ELT(res, 2), 0,
747 COPY_TO_USER_STRING(proj_pj_info(source_crs).definition));
748 SET_STRING_ELT(VECTOR_ELT(res, 3), 0,
749 COPY_TO_USER_STRING(proj_pj_info(target_crs).definition));
750 }
751 SET_STRING_ELT(VECTOR_ELT(res, 4), 0,
752 COPY_TO_USER_STRING(proj_pj_info(pj_transform).definition));
753 }
754
755 /* if ( ob_tran == -1) {
756 for (i=0; i < n; i++) {
757 xx[i] = proj_todeg(xx[i]);
758 yy[i] = proj_todeg(yy[i]);
759 }
760 }*/
761
762 proj_destroy(pj_transform);
763 if (!have_CO) {
764 proj_area_destroy(area_of_interest);
765 proj_destroy(target_crs);
766 proj_destroy(source_crs);
767 }
768 //proj_context_destroy(ctx);
769
770 if (have_z) {
771 for (i=0; i < n; i++) {
772 if (xx[i] == HUGE_VAL || yy[i] == HUGE_VAL || zz[i] == HUGE_VAL
773 || ISNAN(xx[i]) || ISNAN(yy[i]) || ISNAN(zz[i])) {
774 nwarn++;
775 /* Rprintf("transformed point not finite\n");*/
776 }
777 NUMERIC_POINTER(VECTOR_ELT(res, 0))[i] = xx[i];
778 NUMERIC_POINTER(VECTOR_ELT(res, 1))[i] = yy[i];
779 NUMERIC_POINTER(VECTOR_ELT(res, 2))[i] = zz[i];
780 }
781 } else {
782 for (i=0; i < n; i++) {
783 if (xx[i] == HUGE_VAL || yy[i] == HUGE_VAL
784 || ISNAN(xx[i]) || ISNAN(yy[i])) {
785 nwarn++;
786 /* Rprintf("transformed point not finite\n");*/
787 }
788 NUMERIC_POINTER(VECTOR_ELT(res, 0))[i] = xx[i];
789 NUMERIC_POINTER(VECTOR_ELT(res, 1))[i] = yy[i];
790 }
791 }
792
793 if (nwarn > 0) warning("%d projected point(s) not finite", nwarn);
794 UNPROTECT(pc);
795 return(res);
796 }
797
project_ng_coordOp(SEXP proj,SEXP inv,SEXP aoi,SEXP ob_tran)798 SEXP project_ng_coordOp(SEXP proj, SEXP inv, SEXP aoi, SEXP ob_tran) {
799
800 //PJ_CONTEXT *ctx = proj_context_create();
801 PJ *source_crs, *target_crs;
802 PJ* pj_transform = NULL;
803 PJ_AREA *area_of_interest = 0;
804 int use_ob_tran = LOGICAL_POINTER(ob_tran)[0],
805 use_inv, use_aoi=1;
806
807 proj_log_func(PJ_DEFAULT_CTX, NULL, silent_logger);
808
809 /* if (ob_tran == R_NilValue) use_ob_tran = 0;
810 else if (LOGICAL_POINTER(ob_tran)[0] == 1) use_ob_tran = 1;
811 else if (LOGICAL_POINTER(ob_tran)[0] == 0) use_ob_tran = 0;
812 else use_ob_tran = 0;*/
813
814 if (inv == R_NilValue) use_inv = 0;
815 else if (LOGICAL_POINTER(inv)[0] == 1) use_inv = 1;
816 else if (LOGICAL_POINTER(inv)[0] == 0) use_inv = 0;
817 else use_inv = 0;
818
819 if (aoi == R_NilValue) use_aoi = 0;
820
821
822 //Rprintf("target crs input: %s\n", CHAR(STRING_ELT(proj, 0)));
823 if ((target_crs = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(proj, 0)))) == 0) {
824 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
825 //proj_context_destroy(ctx);
826 error("target crs creation failed: %s", errstr);
827 }
828
829 //Rprintf("target crs: %s\n", proj_as_proj_string(PJ_DEFAULT_CTX, target_crs, PJ_PROJ_5, NULL));
830 if (proj_get_type(target_crs) == PJ_TYPE_GEOGRAPHIC_2D_CRS && use_ob_tran) {
831 if ((source_crs = proj_get_source_crs(PJ_DEFAULT_CTX, target_crs)) == 0) {
832 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
833 proj_destroy(target_crs);
834 //proj_context_destroy(ctx);
835 error("source crs creation failed: %s", errstr);
836 }
837 } else {
838 if ((source_crs = proj_crs_get_geodetic_crs(PJ_DEFAULT_CTX, target_crs)) == 0) {
839 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
840 proj_destroy(target_crs);
841 //proj_context_destroy(ctx);
842 error("source crs creation failed: %s", errstr);
843 }
844 }
845 //Rprintf("source crs: %s\n", proj_as_proj_string(PJ_DEFAULT_CTX, source_crs, PJ_PROJ_5, NULL));
846
847 if (use_aoi) {
848 area_of_interest = proj_area_create();
849 proj_area_set_bbox(area_of_interest, NUMERIC_POINTER(aoi)[0],
850 NUMERIC_POINTER(aoi)[1], NUMERIC_POINTER(aoi)[2],
851 NUMERIC_POINTER(aoi)[3]);
852 }
853 //Rprintf("use_inv: %d\n", use_inv);
854 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 2
855 if (use_inv) pj_transform = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
856 proj_as_wkt(PJ_DEFAULT_CTX, target_crs, PJ_WKT2_2018, NULL),
857 proj_as_wkt(PJ_DEFAULT_CTX, source_crs, PJ_WKT2_2018, NULL),
858 area_of_interest);
859 else pj_transform = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
860 proj_as_wkt(PJ_DEFAULT_CTX, source_crs, PJ_WKT2_2018, NULL),
861 proj_as_wkt(PJ_DEFAULT_CTX, target_crs, PJ_WKT2_2018, NULL),
862 area_of_interest);
863 #else
864 if (use_inv) pj_transform = proj_create_crs_to_crs_from_pj(PJ_DEFAULT_CTX, target_crs,
865 source_crs, area_of_interest, NULL);
866 else pj_transform = proj_create_crs_to_crs_from_pj(PJ_DEFAULT_CTX, source_crs,
867 target_crs, area_of_interest, NULL);
868 // FIXME >= 6.2.0
869 #endif
870 if (pj_transform == 0) {
871 proj_area_destroy(area_of_interest);
872 proj_destroy(target_crs);
873 proj_destroy(source_crs);
874 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
875 //proj_context_destroy(ctx);
876 error("coordinate operation creation failed: %s", errstr);
877 }
878 //Rprintf("%s\n", proj_pj_info(pj_transform).definition);
879 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 1
880 warning("no Coordinate Operation normalization available before PROJ 6.1");
881 #else
882 PJ* pj_trans_orig = pj_transform;
883 pj_transform = proj_normalize_for_visualization(PJ_DEFAULT_CTX, pj_trans_orig);
884 proj_destroy(pj_trans_orig);
885 #endif
886 //Rprintf("%s\n", proj_pj_info(pj_transform).definition);
887
888 SEXP res;
889 PROTECT(res = NEW_CHARACTER(1));
890 SET_STRING_ELT(res, 0,
891 COPY_TO_USER_STRING(proj_pj_info(pj_transform).definition));
892 UNPROTECT(1);
893 proj_destroy(pj_transform);
894 proj_area_destroy(area_of_interest);
895 proj_destroy(target_crs);
896 proj_destroy(source_crs);
897 //proj_context_destroy(ctx);
898
899 return(res);
900 }
901
project_ng(SEXP n,SEXP xlon,SEXP ylat,SEXP zz,SEXP coordOp)902 SEXP project_ng(SEXP n, SEXP xlon, SEXP ylat, SEXP zz, SEXP coordOp) {
903 int i, nwarn=0, nn=INTEGER_POINTER(n)[0];
904 SEXP res;
905 //PJ_CONTEXT *ctx = proj_context_create();
906 PJ* pj_transform = NULL;
907 PJ_COORD a, b;
908 double ixlon, iylat, iz=0.0;
909
910 proj_log_func(PJ_DEFAULT_CTX, NULL, silent_logger);
911
912 if ((pj_transform = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(coordOp, 0)))) == 0) {
913 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
914 //proj_context_destroy(ctx);
915 error("coordinate operation creation failed: %s", errstr);
916 }
917 //Rprintf("%s\n", proj_pj_info(pj_transform).definition);
918
919 if (zz != R_NilValue) {
920 PROTECT(res = NEW_LIST(3));
921 SET_VECTOR_ELT(res, 2, NEW_NUMERIC(nn));
922 } else {
923 PROTECT(res = NEW_LIST(2));
924 }
925 SET_VECTOR_ELT(res, 0, NEW_NUMERIC(nn));
926 SET_VECTOR_ELT(res, 1, NEW_NUMERIC(nn));
927
928 for (i=0; i<nn; i++) {
929 ixlon = NUMERIC_POINTER(xlon)[i];
930 iylat = NUMERIC_POINTER(ylat)[i];
931 if (zz != R_NilValue) iz = NUMERIC_POINTER(zz)[i];
932 // preserve NAs and NaNs. Allow Infs, since maybe proj can handle them.
933 if (ISNAN(ixlon) || ISNAN(iylat)) {
934 NUMERIC_POINTER(VECTOR_ELT(res, 0))[i]=ixlon;
935 NUMERIC_POINTER(VECTOR_ELT(res, 1))[i]=iylat;
936 } else {
937 a = proj_coord(ixlon, iylat, iz, 0);
938 b = proj_trans(pj_transform, PJ_FWD, a);
939 if (b.uv.u == HUGE_VAL || ISNAN(b.uv.u) || b.uv.v == HUGE_VAL ||
940 ISNAN(b.uv.v)) {
941 nwarn++;
942 }
943 NUMERIC_POINTER(VECTOR_ELT(res, 0))[i]=b.uv.u;
944 NUMERIC_POINTER(VECTOR_ELT(res, 1))[i]=b.uv.v;
945 if (zz != R_NilValue)
946 NUMERIC_POINTER(VECTOR_ELT(res, 2))[i]=b.uvw.w;
947 }
948 }
949 if (nwarn > 0) warning("%d projected point(s) not finite", nwarn);
950
951 proj_destroy(pj_transform);
952 //proj_context_destroy(ctx);
953
954 UNPROTECT(1);
955 return(res);
956 }
957
P6_SRID_proj(SEXP inSRID,SEXP format,SEXP multiline,SEXP in_format,SEXP epsg,SEXP out_format)958 SEXP P6_SRID_proj(SEXP inSRID, SEXP format, SEXP multiline, SEXP in_format,
959 SEXP epsg, SEXP out_format) {
960
961 SEXP ans;
962 SEXP Datum, Ellps;
963 int pc=0;
964 int vis_order;
965 SEXP enforce_xy = getAttrib(in_format, install("enforce_xy"));
966 const char *pszSRS = NULL;
967
968 if (enforce_xy == R_NilValue) vis_order = 0;
969 else if (LOGICAL_POINTER(enforce_xy)[0] == 1) vis_order = 1;
970 else if (LOGICAL_POINTER(enforce_xy)[0] == 0) vis_order = 0;
971 else vis_order = 0;
972
973 PJ *source_crs = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(inSRID, 0)));
974 // valgrind resolved by copying EJP
975 // https://github.com/r-spatial/gstat/issues/82#issuecomment-762970800
976 if (source_crs == NULL) {
977 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
978 // proj_context_destroy(ctx);
979 error("source crs creation failed: %s", errstr);
980 }
981
982 #if (PROJ_VERSION_MAJOR > 6 || (PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR >= 1))
983 PJ_TYPE type;
984 type = proj_get_type(source_crs);
985
986 if (type == PJ_TYPE_BOUND_CRS) {
987 PJ *orig_crs = source_crs;
988 source_crs = proj_get_source_crs(PJ_DEFAULT_CTX, orig_crs);
989 proj_destroy(orig_crs);
990 if (source_crs == NULL) {
991 // proj_context_destroy(ctx);
992 error("crs not converted to source only");
993 }
994 }
995 #endif
996
997 #if (PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 3)
998 warning("no CRS normalization available before PROJ 6.3");
999 #else
1000 if (vis_order) {
1001 PJ *orig_crs = source_crs;
1002 source_crs = proj_normalize_for_visualization(PJ_DEFAULT_CTX, orig_crs);
1003 proj_destroy(orig_crs);
1004 if (source_crs == NULL) {
1005 // proj_context_destroy(ctx);
1006 error("crs not converted to visualization order");
1007 }
1008 }
1009 #endif
1010 // FIXME PROJ 8.0 datum ensemble vulnerability
1011 PJ* dtm = proj_crs_get_datum(PJ_DEFAULT_CTX, source_crs);
1012 if (dtm != NULL) {
1013 PROTECT(Datum = NEW_CHARACTER(1)); pc++;
1014 SET_STRING_ELT(Datum, 0, COPY_TO_USER_STRING(proj_get_name(dtm)));
1015 proj_destroy(dtm);
1016 } else {
1017 Datum = R_NilValue;
1018 }
1019
1020 PJ* ellps = proj_get_ellipsoid(PJ_DEFAULT_CTX, source_crs);
1021 if (ellps != NULL) {
1022 PROTECT(Ellps = NEW_CHARACTER(1)); pc++;
1023 SET_STRING_ELT(Ellps, 0, COPY_TO_USER_STRING(proj_get_name(ellps)));
1024 proj_destroy(ellps);
1025 } else {
1026 Ellps = R_NilValue;
1027 }
1028
1029
1030
1031 PROTECT(ans=NEW_CHARACTER(1)); pc++;
1032
1033 if (INTEGER_POINTER(out_format)[0] == 1L) {
1034
1035 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 3
1036 if ((pszSRS = proj_as_wkt(PJ_DEFAULT_CTX, source_crs, PJ_WKT2_2018, NULL))
1037 == NULL) {
1038 #else
1039 if ((pszSRS = proj_as_wkt(PJ_DEFAULT_CTX, source_crs, PJ_WKT2_2019, NULL))
1040 == NULL) {
1041 #endif
1042 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
1043 warning("export to WKT2 failed: %s", errstr);
1044 SET_STRING_ELT(ans, 0, NA_STRING);
1045 } else {
1046 SET_STRING_ELT(ans, 0, COPY_TO_USER_STRING(pszSRS));
1047 }
1048 } else if (INTEGER_POINTER(out_format)[0] == 2L) {
1049
1050 if ((pszSRS = proj_as_proj_string(PJ_DEFAULT_CTX, source_crs, PJ_PROJ_4, NULL))
1051 == NULL) {
1052 const char *errstr = proj_errno_string(proj_context_errno(PJ_DEFAULT_CTX));
1053 warning("export to PROJ failed: %s", errstr);
1054 SET_STRING_ELT(ans, 0, NA_STRING);
1055 } else {
1056 SET_STRING_ELT(ans, 0, COPY_TO_USER_STRING(pszSRS));
1057 }
1058 } else {
1059 proj_destroy(source_crs);
1060 // proj_context_destroy(ctx);
1061 UNPROTECT(pc);
1062 error("unknown output format");
1063 }
1064 setAttrib(ans, install("datum"), Datum);
1065 setAttrib(ans, install("ellps"), Ellps);
1066
1067
1068 proj_destroy(source_crs);
1069 // proj_context_destroy(ctx);
1070 UNPROTECT(pc);
1071
1072 return(ans);
1073
1074 }
1075
1076
1077
1078 // blocks error messages in this context
1079 // https://lists.osgeo.org/pipermail/proj/2019-March/008310.html
1080 static void proj_logger(void * user_data, int level, const char * message) {}
1081
1082 // code borrowed from GRASS g.proj main.c adapted for PROJ6 by Markus Metz
1083
1084 SEXP
1085 PROJcopyEPSG(SEXP tf) {
1086
1087 SEXP ans;
1088 PROTECT(ans=NEW_INTEGER(1));
1089 INTEGER_POINTER(ans)[0] = 0;
1090 int i, crs_cnt;
1091 PROJ_CRS_INFO **proj_crs_info;
1092 PJ_CONTEXT *ctx = proj_context_create();
1093 FILE *fptf;
1094
1095
1096 crs_cnt = 0;
1097 proj_crs_info = proj_get_crs_info_list_from_database(ctx, "EPSG", NULL,
1098 &crs_cnt);
1099 if (crs_cnt < 1) {
1100 UNPROTECT(1);
1101 return(ans);
1102 }
1103 fptf = fopen(CHAR(STRING_ELT(tf, 0)), "wb");
1104 if (fptf == NULL) {
1105 UNPROTECT(1);
1106 return(ans);
1107 }
1108 fprintf(fptf, "\"code\",\"note\",\"prj4\",\"prj_method\"\n");
1109
1110 // PJ *pj = NULL;
1111 // blocks error messages in this context
1112 proj_log_func(ctx, NULL, proj_logger);
1113 for (i = 0; i < crs_cnt; i++) {
1114 const char *proj_definition;
1115
1116 PJ* pj = proj_create_from_database(ctx,
1117 proj_crs_info[i]->auth_name, proj_crs_info[i]->code,
1118 PJ_CATEGORY_CRS, 0, NULL);
1119 proj_definition = proj_as_proj_string(ctx, pj, // valgrind
1120 PJ_PROJ_5, NULL);
1121 // proj_destroy(pj); // valgrind
1122
1123 fprintf(fptf, "%s,\"%s\",\"%s\",\"%s\"\n", proj_crs_info[i]->code, //ASAN
1124 proj_crs_info[i]->name, proj_definition,
1125 proj_crs_info[i]->projection_method_name);
1126 proj_destroy(pj);
1127 }
1128
1129 fclose(fptf);
1130 proj_crs_info_list_destroy(proj_crs_info);
1131 proj_context_destroy(ctx);
1132 INTEGER_POINTER(ans)[0] = crs_cnt;
1133 UNPROTECT(1);
1134
1135 return(ans);
1136
1137 }
1138
1139 #ifdef __cplusplus
1140 }
1141 #endif
1142 #endif // PROJ_H_API
1143
1144
1145