1 #include <grass/gis.h>
2 #include <grass/gprojects.h>
3 #include <grass/glocale.h>
4
5 #include <ogr_srs_api.h>
6 #include <cpl_conv.h>
7 #include "global.h"
8
9 /* get projection info of OGR layer in GRASS format
10 * return 0 on success (some non-xy SRS)
11 * return 1 if no SRS available
12 * return 2 if SRS available but unreadable */
get_layer_proj(OGRLayerH Ogr_layer,struct Cell_head * cellhd,struct Key_Value ** proj_info,struct Key_Value ** proj_units,char ** proj_srid,char ** proj_wkt,char * geom_col,int verbose)13 int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
14 struct Key_Value **proj_info,
15 struct Key_Value **proj_units,
16 char **proj_srid, char **proj_wkt,
17 char *geom_col, int verbose)
18 {
19 OGRSpatialReferenceH hSRS;
20
21 hSRS = NULL;
22 *proj_info = NULL;
23 *proj_units = NULL;
24 *proj_srid = NULL;
25 *proj_wkt = NULL;
26
27 /* Fetch input layer projection in GRASS form. */
28 #if GDAL_VERSION_NUM >= 1110000
29 if (geom_col) {
30 int igeom;
31 OGRGeomFieldDefnH Ogr_geomdefn;
32 OGRFeatureDefnH Ogr_featuredefn;
33
34 Ogr_featuredefn = OGR_L_GetLayerDefn(Ogr_layer);
35 igeom = OGR_FD_GetGeomFieldIndex(Ogr_featuredefn, geom_col);
36 if (igeom < 0)
37 G_fatal_error(_("Geometry column <%s> not found in input layer <%s>"),
38 geom_col, OGR_L_GetName(Ogr_layer));
39 Ogr_geomdefn = OGR_FD_GetGeomFieldDefn(Ogr_featuredefn, igeom);
40 hSRS = OGR_GFld_GetSpatialRef(Ogr_geomdefn);
41 }
42 else {
43 hSRS = OGR_L_GetSpatialRef(Ogr_layer);
44 }
45 #else
46 hSRS = OGR_L_GetSpatialRef(Ogr_layer); /* should not be freed later */
47 #endif
48
49 /* verbose is used only when comparing input SRS to GRASS projection,
50 * not when comparing SRS's of several input layers */
51 if (GPJ_osr_to_grass(cellhd, proj_info,
52 proj_units, hSRS, 0) < 0) {
53 /* TODO: GPJ_osr_to_grass() does not return anything < 0
54 * check with GRASS 6 and GRASS 5 */
55 G_warning(_("Unable to convert input layer projection information to "
56 "GRASS format for checking"));
57 if (verbose && hSRS != NULL) {
58 char *wkt = NULL;
59
60 if (OSRExportToPrettyWkt(hSRS, &wkt, FALSE) != OGRERR_NONE) {
61 G_warning(_("Can't get WKT parameter string"));
62 }
63 else if (wkt) {
64 G_important_message(_("WKT definition:\n%s"), wkt);
65 }
66 }
67
68 return 2;
69 }
70 /* custom checks because if in doubt GPJ_osr_to_grass() returns a
71 * xy CRS */
72 if (hSRS == NULL) {
73 if (verbose) {
74 G_important_message(_("No projection information available for layer <%s>"),
75 OGR_L_GetName(Ogr_layer));
76 }
77
78 return 1;
79 }
80
81 if (!OSRIsProjected(hSRS) && !OSRIsGeographic(hSRS)) {
82 G_important_message(_("Projection for layer <%s> does not contain a valid SRS"),
83 OGR_L_GetName(Ogr_layer));
84
85 if (verbose) {
86 char *wkt = NULL;
87
88 if (OSRExportToPrettyWkt(hSRS, &wkt, FALSE) != OGRERR_NONE) {
89 G_important_message(_("Can't get WKT parameter string"));
90 }
91 else if (wkt) {
92 G_important_message(_("WKT definition:\n%s"), wkt);
93 }
94 }
95
96 return 2;
97 }
98 else{
99 const char *authkey, *authname, *authcode;
100 #if GDAL_VERSION_NUM >= 3000000
101 char **papszOptions;
102
103 /* get WKT2 definition */
104 papszOptions = G_calloc(3, sizeof(char *));
105 papszOptions[0] = G_store("MULTILINE=YES");
106 papszOptions[1] = G_store("FORMAT=WKT2");
107 OSRExportToWktEx(hSRS, proj_wkt, (const char **)papszOptions);
108 G_free(papszOptions[0]);
109 G_free(papszOptions[1]);
110 G_free(papszOptions);
111 #else
112 OSRExportToWkt(hSRS, proj_wkt);
113 #endif
114
115 if (OSRIsProjected(hSRS))
116 authkey = "PROJCS";
117 else /* is geographic */
118 authkey = "GEOGCS";
119
120 authname = OSRGetAuthorityName(hSRS, authkey);
121 if (authname && *authname) {
122 authcode = OSRGetAuthorityCode(hSRS, authkey);
123 if (authcode && *authcode) {
124 G_asprintf(proj_srid, "%s:%s", authname, authcode);
125 }
126 }
127 }
128
129 return 0;
130 }
131
132
133 /* compare projections of all OGR layers
134 * return 0 if all layers have the same projection
135 * return 1 if layer projections differ */
cmp_layer_srs(ds_t Ogr_ds,int nlayers,int * layers,char ** layer_names,char * geom_col)136 int cmp_layer_srs(ds_t Ogr_ds, int nlayers, int *layers,
137 char **layer_names, char *geom_col)
138 {
139 int layer;
140 struct Key_Value *proj_info1, *proj_units1;
141 char *proj_srid1, *proj_wkt1;
142 struct Key_Value *proj_info2, *proj_units2;
143 char *proj_srid2, *proj_wkt2;
144 struct Cell_head cellhd1, cellhd2;
145 OGRLayerH Ogr_layer;
146
147 if (nlayers == 1)
148 return 0;
149
150 proj_info1 = proj_units1 = NULL;
151 proj_srid1 = proj_wkt1 = NULL;
152 proj_info2 = proj_units2 = NULL;
153 proj_srid2 = proj_wkt2 = NULL;
154
155 G_get_window(&cellhd1);
156 layer = 0;
157 do {
158 /* Get first SRS */
159 Ogr_layer = ds_getlayerbyindex(Ogr_ds, layers[layer]);
160
161 if (get_layer_proj(Ogr_layer, &cellhd1, &proj_info1, &proj_units1,
162 &proj_srid1, &proj_wkt1, geom_col, 0) == 0) {
163 break;
164 }
165 layer++;
166 } while (layer < nlayers);
167
168 if (layer == nlayers) {
169 /* could not get layer proj in GRASS format for any of the layers
170 * -> projections of all layers are the same, i.e. unreadable by GRASS */
171 G_warning(_("Layer projections are unreadable"));
172 if (proj_info1)
173 G_free_key_value(proj_info1);
174 if (proj_units1)
175 G_free_key_value(proj_units1);
176 if (proj_srid1)
177 G_free(proj_srid1);
178 if (proj_wkt1)
179 CPLFree(proj_wkt1);
180
181 return 0;
182 }
183 if (layer > 0) {
184 /* could not get layer proj in GRASS format for at least one of the layers
185 * -> mix of unreadable and readable projections */
186 G_warning(_("Projection for layer <%s> is unreadable"),
187 layer_names[layer]);
188 if (proj_info1)
189 G_free_key_value(proj_info1);
190 if (proj_units1)
191 G_free_key_value(proj_units1);
192 if (proj_srid1)
193 G_free(proj_srid1);
194 if (proj_wkt1)
195 CPLFree(proj_wkt1);
196
197 return 1;
198 }
199
200 for (layer = 1; layer < nlayers; layer++) {
201 /* Get SRS of other layer(s) */
202 Ogr_layer = ds_getlayerbyindex(Ogr_ds, layers[layer]);
203 G_get_window(&cellhd2);
204 if (get_layer_proj(Ogr_layer, &cellhd2, &proj_info2, &proj_units2,
205 &proj_srid2, &proj_wkt2, geom_col, 0) != 0) {
206 if (proj_info1)
207 G_free_key_value(proj_info1);
208 if (proj_units1)
209 G_free_key_value(proj_units1);
210 if (proj_srid1)
211 G_free(proj_srid1);
212 if (proj_wkt1)
213 CPLFree(proj_wkt1);
214
215 return 1;
216 }
217
218 if (cellhd1.proj != cellhd2.proj
219 || G_compare_projections(proj_info1, proj_units1,
220 proj_info2, proj_units2) < 0) {
221 if (proj_info1)
222 G_free_key_value(proj_info1);
223 if (proj_units1)
224 G_free_key_value(proj_units1);
225 if (proj_srid1)
226 G_free(proj_srid1);
227 if (proj_wkt1)
228 CPLFree(proj_wkt1);
229 if (proj_info2)
230 G_free_key_value(proj_info2);
231 if (proj_units2)
232 G_free_key_value(proj_units2);
233 if (proj_srid2)
234 G_free(proj_srid2);
235 if (proj_wkt2)
236 CPLFree(proj_wkt2);
237
238 G_warning(_("Projection of layer <%s> is different from "
239 "projection of layer <%s>"),
240 layer_names[layer], layer_names[layer - 1]);
241
242 return 1;
243 }
244 if (proj_info2)
245 G_free_key_value(proj_info2);
246 if (proj_units2)
247 G_free_key_value(proj_units2);
248 if (proj_srid2)
249 G_free(proj_srid2);
250 if (proj_wkt2)
251 CPLFree(proj_wkt2);
252 }
253 if (proj_info1)
254 G_free_key_value(proj_info1);
255 if (proj_units1)
256 G_free_key_value(proj_units1);
257 if (proj_srid1)
258 G_free(proj_srid1);
259 if (proj_wkt1)
260 CPLFree(proj_wkt1);
261
262 return 0;
263 }
264
265 /* keep in sync with r.in.gdal, r.external, v.external */
check_projection(struct Cell_head * cellhd,ds_t hDS,int layer,char * geom_col,char * outloc,int create_only,int override,int check_only)266 void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_col,
267 char *outloc, int create_only, int override,
268 int check_only)
269 {
270 struct Cell_head loc_wind;
271 struct Key_Value *proj_info = NULL,
272 *proj_units = NULL;
273 struct Key_Value *loc_proj_info = NULL,
274 *loc_proj_units = NULL;
275 char *wkt = NULL, *srid = NULL;
276 char error_msg[8096];
277 int proj_trouble;
278 OGRLayerH Ogr_layer;
279
280 /* Get first layer to be imported to use for projection check */
281 Ogr_layer = ds_getlayerbyindex(hDS, layer);
282
283 /* -------------------------------------------------------------------- */
284 /* Fetch the projection in GRASS form, SRID, and WKT. */
285 /* -------------------------------------------------------------------- */
286
287 /* proj_trouble:
288 * 0: valid srs
289 * 1: no srs, default to xy
290 * 2: unreadable srs, default to xy
291 */
292
293 /* Projection only required for checking so convert non-interactively */
294 proj_trouble = get_layer_proj(Ogr_layer, cellhd, &proj_info, &proj_units,
295 &srid, &wkt, geom_col, 1);
296
297 /* -------------------------------------------------------------------- */
298 /* Do we need to create a new location? */
299 /* -------------------------------------------------------------------- */
300 if (outloc != NULL) {
301 /* do not create a xy location because this can mean that the
302 * real SRS has not been recognized or is missing */
303 if (proj_trouble) {
304 G_fatal_error(_("Unable to convert input map projection to GRASS "
305 "format; cannot create new location."));
306 }
307 else {
308 if (0 != G_make_location_crs(outloc, cellhd, proj_info,
309 proj_units, srid, wkt)) {
310 G_fatal_error(_("Unable to create new location <%s>"),
311 outloc);
312 }
313 G_message(_("Location <%s> created"), outloc);
314
315 G_unset_window(); /* new location, projection, and window */
316 G_get_window(cellhd);
317 }
318
319 /* If create only, clean up and exit here */
320 if (create_only) {
321 ds_close(hDS);
322 exit(EXIT_SUCCESS);
323 }
324 }
325 else {
326 int err = 0;
327 void (*msg_fn)(const char *, ...);
328
329 if (check_only && override) {
330 /* can't check when over-riding check */
331 override = 0;
332 }
333
334 if (proj_trouble == 2) {
335 strcpy(error_msg, _("Unable to convert input map projection information "
336 "to GRASS format."));
337 if (override) {
338 msg_fn = G_warning;
339 }
340 else {
341 msg_fn = G_fatal_error;
342 ds_close(hDS);
343 }
344 msg_fn(error_msg);
345 if (!override) {
346 exit(EXIT_FAILURE);
347 }
348 }
349
350 /* -------------------------------------------------------------------- */
351 /* Does the projection of the current location match the */
352 /* dataset? */
353 /* -------------------------------------------------------------------- */
354 G_get_default_window(&loc_wind);
355 /* fetch LOCATION PROJ info */
356 if (loc_wind.proj != PROJECTION_XY) {
357 loc_proj_info = G_get_projinfo();
358 loc_proj_units = G_get_projunits();
359 }
360
361 if (override) {
362 cellhd->proj = loc_wind.proj;
363 cellhd->zone = loc_wind.zone;
364 G_message(_("Over-riding projection check"));
365 }
366 else if (loc_wind.proj != cellhd->proj
367 || (err =
368 G_compare_projections(loc_proj_info, loc_proj_units,
369 proj_info, proj_units)) != 1) {
370 int i_value;
371
372 strcpy(error_msg,
373 _("Projection of dataset does not"
374 " appear to match current location.\n\n"));
375
376 /* TODO: output this info sorted by key: */
377 if (loc_wind.proj != cellhd->proj || err != -2) {
378 /* error in proj_info */
379 if (loc_proj_info != NULL) {
380 strcat(error_msg, _("Location PROJ_INFO is:\n"));
381 for (i_value = 0; i_value < loc_proj_info->nitems;
382 i_value++)
383 sprintf(error_msg + strlen(error_msg), "%s: %s\n",
384 loc_proj_info->key[i_value],
385 loc_proj_info->value[i_value]);
386 strcat(error_msg, "\n");
387 }
388 else {
389 strcat(error_msg, _("Location PROJ_INFO is:\n"));
390 if (loc_wind.proj == PROJECTION_XY)
391 sprintf(error_msg + strlen(error_msg),
392 "Location proj = %d (unreferenced/unknown)\n",
393 loc_wind.proj);
394 else if (loc_wind.proj == PROJECTION_LL)
395 sprintf(error_msg + strlen(error_msg),
396 "Location proj = %d (lat/long)\n",
397 loc_wind.proj);
398 else if (loc_wind.proj == PROJECTION_UTM)
399 sprintf(error_msg + strlen(error_msg),
400 "Location proj = %d (UTM), zone = %d\n",
401 loc_wind.proj, cellhd->zone);
402 else
403 sprintf(error_msg + strlen(error_msg),
404 "Location proj = %d (unknown), zone = %d\n",
405 loc_wind.proj, cellhd->zone);
406 }
407
408 if (proj_info != NULL) {
409 strcat(error_msg, _("Dataset PROJ_INFO is:\n"));
410 for (i_value = 0; i_value < proj_info->nitems; i_value++)
411 sprintf(error_msg + strlen(error_msg), "%s: %s\n",
412 proj_info->key[i_value],
413 proj_info->value[i_value]);
414 }
415 else {
416 strcat(error_msg, _("Dataset PROJ_INFO is:\n"));
417 if (cellhd->proj == PROJECTION_XY)
418 sprintf(error_msg + strlen(error_msg),
419 "Dataset proj = %d (unreferenced/unknown)\n",
420 cellhd->proj);
421 else if (cellhd->proj == PROJECTION_LL)
422 sprintf(error_msg + strlen(error_msg),
423 "Dataset proj = %d (lat/long)\n",
424 cellhd->proj);
425 else if (cellhd->proj == PROJECTION_UTM)
426 sprintf(error_msg + strlen(error_msg),
427 "Dataset proj = %d (UTM), zone = %d\n",
428 cellhd->proj, cellhd->zone);
429 else
430 sprintf(error_msg + strlen(error_msg),
431 "Dataset proj = %d (unknown), zone = %d\n",
432 cellhd->proj, cellhd->zone);
433 }
434 if (loc_wind.proj != cellhd->proj) {
435 strcat(error_msg, "\nDifference in: proj\n");
436 }
437 else {
438 strcat(error_msg, "\nDifference in: ");
439 switch (err) {
440 case -1:
441 strcat(error_msg, "proj\n");
442 break;
443 case -2:
444 strcat(error_msg, "units\n");
445 break;
446 case -3:
447 strcat(error_msg, "datum\n");
448 break;
449 case -4:
450 strcat(error_msg, "ellps, a, es\n");
451 break;
452 case -5:
453 strcat(error_msg, "zone\n");
454 break;
455 case -6:
456 strcat(error_msg, "south\n");
457 break;
458 case -7:
459 strcat(error_msg, "x_0\n");
460 break;
461 case -8:
462 strcat(error_msg, "y_0\n");
463 break;
464 case -9:
465 strcat(error_msg, "lon_0\n");
466 break;
467 case -10:
468 strcat(error_msg, "lat_0\n");
469 break;
470 case -11:
471 strcat(error_msg, "lat_1, lat2\n");
472 break;
473 }
474 }
475 }
476 else {
477 /* error in proj_units */
478 if (loc_proj_units != NULL) {
479 strcat(error_msg, "Location PROJ_UNITS is:\n");
480 for (i_value = 0; i_value < loc_proj_units->nitems;
481 i_value++)
482 sprintf(error_msg + strlen(error_msg), "%s: %s\n",
483 loc_proj_units->key[i_value],
484 loc_proj_units->value[i_value]);
485 strcat(error_msg, "\n");
486 }
487
488 if (proj_units != NULL) {
489 strcat(error_msg, "Dataset PROJ_UNITS is:\n");
490 for (i_value = 0; i_value < proj_units->nitems; i_value++)
491 sprintf(error_msg + strlen(error_msg), "%s: %s\n",
492 proj_units->key[i_value],
493 proj_units->value[i_value]);
494 }
495 }
496 if (!check_only) {
497 strcat(error_msg,
498 _("\nIn case of no significant differences in the projection definitions,"
499 " use the -o flag to ignore them and use"
500 " current location definition.\n"));
501 strcat(error_msg,
502 _("Consider generating a new location from the input dataset using "
503 "the 'location' parameter.\n"));
504 }
505
506 if (check_only)
507 msg_fn = G_message;
508 else
509 msg_fn = G_fatal_error;
510 msg_fn("%s", error_msg);
511 if (check_only) {
512 ds_close(hDS);
513 exit(EXIT_FAILURE);
514 }
515 }
516 else {
517 if (check_only)
518 msg_fn = G_message;
519 else
520 msg_fn = G_verbose_message;
521 msg_fn(_("Projection of input dataset and current location "
522 "appear to match"));
523 if (check_only) {
524 ds_close(hDS);
525 exit(EXIT_SUCCESS);
526 }
527 }
528 }
529 }
530