1 #include <grass/gis.h>
2 #include <grass/gprojects.h>
3 #include <grass/glocale.h>
4
5 #include <ogr_srs_api.h>
6 #include "local_proto.h"
7
8 /* get projection info of OGR layer in GRASS format
9 * return 0 on success (some non-xy SRS)
10 * return 1 if no SRS available
11 * 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)12 int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
13 struct Key_Value **proj_info,
14 struct Key_Value **proj_units,
15 char **proj_srid, char **proj_wkt,
16 char *geom_col, int verbose)
17 {
18 OGRSpatialReferenceH hSRS;
19
20 hSRS = NULL;
21 *proj_info = NULL;
22 *proj_units = NULL;
23 *proj_srid = NULL;
24 *proj_wkt = NULL;
25
26 /* Fetch input layer projection in GRASS form. */
27 #if GDAL_VERSION_NUM >= 1110000
28 if (geom_col) {
29 int igeom;
30 OGRGeomFieldDefnH Ogr_geomdefn;
31 OGRFeatureDefnH Ogr_featuredefn;
32
33 Ogr_featuredefn = OGR_L_GetLayerDefn(Ogr_layer);
34 igeom = OGR_FD_GetGeomFieldIndex(Ogr_featuredefn, geom_col);
35 if (igeom < 0)
36 G_fatal_error(_("Geometry column <%s> not found in input layer <%s>"),
37 geom_col, OGR_L_GetName(Ogr_layer));
38 Ogr_geomdefn = OGR_FD_GetGeomFieldDefn(Ogr_featuredefn, igeom);
39 hSRS = OGR_GFld_GetSpatialRef(Ogr_geomdefn);
40 }
41 else {
42 hSRS = OGR_L_GetSpatialRef(Ogr_layer);
43 }
44 #else
45 hSRS = OGR_L_GetSpatialRef(Ogr_layer); /* should not be freed later */
46 #endif
47
48 /* verbose is used only when comparing input SRS to GRASS projection,
49 * not when comparing SRS's of several input layers */
50 if (GPJ_osr_to_grass(cellhd, proj_info,
51 proj_units, hSRS, 0) < 0) {
52 /* TODO: GPJ_osr_to_grass() does not return anything < 0
53 * check with GRASS 6 and GRASS 5 */
54 G_warning(_("Unable to convert input layer projection information to "
55 "GRASS format for checking"));
56 if (verbose && hSRS != NULL) {
57 char *wkt = NULL;
58
59 if (OSRExportToPrettyWkt(hSRS, &wkt, FALSE) != OGRERR_NONE) {
60 G_warning(_("Can't get WKT parameter string"));
61 }
62 else if (wkt) {
63 G_important_message(_("WKT definition:\n%s"), wkt);
64 }
65 }
66
67 return 2;
68 }
69 /* custom checks because if in doubt GPJ_osr_to_grass() returns a
70 * xy CRS */
71 if (hSRS == NULL) {
72 if (verbose) {
73 G_important_message(_("No projection information available for layer <%s>"),
74 OGR_L_GetName(Ogr_layer));
75 }
76
77 return 1;
78 }
79
80 if (!OSRIsProjected(hSRS) && !OSRIsGeographic(hSRS)) {
81 G_important_message(_("Projection for layer <%s> does not contain a valid SRS"),
82 OGR_L_GetName(Ogr_layer));
83
84 if (verbose) {
85 char *wkt = NULL;
86
87 if (OSRExportToPrettyWkt(hSRS, &wkt, FALSE) != OGRERR_NONE) {
88 G_important_message(_("Can't get WKT parameter string"));
89 }
90 else if (wkt) {
91 G_important_message(_("WKT definition:\n%s"), wkt);
92 }
93 }
94
95 return 2;
96 }
97 else{
98 const char *authkey, *authname, *authcode;
99 #if GDAL_VERSION_NUM >= 3000000
100 char **papszOptions;
101
102 /* get WKT2 definition */
103 papszOptions = G_calloc(3, sizeof(char *));
104 papszOptions[0] = G_store("MULTILINE=YES");
105 papszOptions[1] = G_store("FORMAT=WKT2");
106 OSRExportToWktEx(hSRS, proj_wkt, (const char **)papszOptions);
107 G_free(papszOptions[0]);
108 G_free(papszOptions[1]);
109 G_free(papszOptions);
110 #else
111 OSRExportToWkt(hSRS, proj_wkt);
112 #endif
113
114 if (OSRIsProjected(hSRS))
115 authkey = "PROJCS";
116 else /* is geographic */
117 authkey = "GEOGCS";
118
119 authname = OSRGetAuthorityName(hSRS, authkey);
120 if (authname && *authname) {
121 authcode = OSRGetAuthorityCode(hSRS, authkey);
122 if (authcode && *authcode) {
123 G_asprintf(proj_srid, "%s:%s", authname, authcode);
124 }
125 }
126 }
127
128 return 0;
129 }
130
131 /* keep in sync with r.in.gdal, r.external, v.in.ogr */
check_projection(struct Cell_head * cellhd,ds_t hDS,int layer,char * geom_col,char * outloc,int create_only,int override,int check_only)132 void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_col,
133 char *outloc, int create_only, int override,
134 int check_only)
135 {
136 struct Cell_head loc_wind;
137 struct Key_Value *proj_info = NULL,
138 *proj_units = NULL;
139 struct Key_Value *loc_proj_info = NULL,
140 *loc_proj_units = NULL;
141 char *wkt = NULL, *srid = NULL;
142 char error_msg[8096];
143 int proj_trouble;
144 OGRLayerH Ogr_layer;
145
146 /* Get first layer to be imported to use for projection check */
147 Ogr_layer = ds_getlayerbyindex(hDS, layer);
148
149 /* -------------------------------------------------------------------- */
150 /* Fetch the projection in GRASS form, SRID, and WKT. */
151 /* -------------------------------------------------------------------- */
152
153 /* proj_trouble:
154 * 0: valid srs
155 * 1: no srs, default to xy
156 * 2: unreadable srs, default to xy
157 */
158
159 /* Projection only required for checking so convert non-interactively */
160 proj_trouble = get_layer_proj(Ogr_layer, cellhd, &proj_info, &proj_units,
161 &srid, &wkt, geom_col, 1);
162
163 /* -------------------------------------------------------------------- */
164 /* Do we need to create a new location? */
165 /* -------------------------------------------------------------------- */
166 if (outloc != NULL) {
167 /* do not create a xy location because this can mean that the
168 * real SRS has not been recognized or is missing */
169 if (proj_trouble) {
170 G_fatal_error(_("Unable to convert input map projection to GRASS "
171 "format; cannot create new location."));
172 }
173 else {
174 if (0 != G_make_location_crs(outloc, cellhd, proj_info,
175 proj_units, srid, wkt)) {
176 G_fatal_error(_("Unable to create new location <%s>"),
177 outloc);
178 }
179 G_message(_("Location <%s> created"), outloc);
180
181 G_unset_window(); /* new location, projection, and window */
182 G_get_window(cellhd);
183 }
184
185 /* If create only, clean up and exit here */
186 if (create_only) {
187 ds_close(hDS);
188 exit(EXIT_SUCCESS);
189 }
190 }
191 else {
192 int err = 0;
193 void (*msg_fn)(const char *, ...);
194
195 if (check_only && override) {
196 /* can't check when over-riding check */
197 override = 0;
198 }
199
200 if (proj_trouble == 2) {
201 strcpy(error_msg, _("Unable to convert input map projection information "
202 "to GRASS format."));
203 if (override) {
204 msg_fn = G_warning;
205 }
206 else {
207 msg_fn = G_fatal_error;
208 ds_close(hDS);
209 }
210 msg_fn(error_msg);
211 if (!override) {
212 exit(EXIT_FAILURE);
213 }
214 }
215
216 /* -------------------------------------------------------------------- */
217 /* Does the projection of the current location match the */
218 /* dataset? */
219 /* -------------------------------------------------------------------- */
220 G_get_default_window(&loc_wind);
221 /* fetch LOCATION PROJ info */
222 if (loc_wind.proj != PROJECTION_XY) {
223 loc_proj_info = G_get_projinfo();
224 loc_proj_units = G_get_projunits();
225 }
226
227 if (override) {
228 cellhd->proj = loc_wind.proj;
229 cellhd->zone = loc_wind.zone;
230 G_message(_("Over-riding projection check"));
231 }
232 else if (loc_wind.proj != cellhd->proj
233 || (err =
234 G_compare_projections(loc_proj_info, loc_proj_units,
235 proj_info, proj_units)) != 1) {
236 int i_value;
237
238 strcpy(error_msg,
239 _("Projection of dataset does not"
240 " appear to match current location.\n\n"));
241
242 /* TODO: output this info sorted by key: */
243 if (loc_wind.proj != cellhd->proj || err != -2) {
244 /* error in proj_info */
245 if (loc_proj_info != NULL) {
246 strcat(error_msg, _("Location PROJ_INFO is:\n"));
247 for (i_value = 0; i_value < loc_proj_info->nitems;
248 i_value++)
249 sprintf(error_msg + strlen(error_msg), "%s: %s\n",
250 loc_proj_info->key[i_value],
251 loc_proj_info->value[i_value]);
252 strcat(error_msg, "\n");
253 }
254 else {
255 strcat(error_msg, _("Location PROJ_INFO is:\n"));
256 if (loc_wind.proj == PROJECTION_XY)
257 sprintf(error_msg + strlen(error_msg),
258 "Location proj = %d (unreferenced/unknown)\n",
259 loc_wind.proj);
260 else if (loc_wind.proj == PROJECTION_LL)
261 sprintf(error_msg + strlen(error_msg),
262 "Location proj = %d (lat/long)\n",
263 loc_wind.proj);
264 else if (loc_wind.proj == PROJECTION_UTM)
265 sprintf(error_msg + strlen(error_msg),
266 "Location proj = %d (UTM), zone = %d\n",
267 loc_wind.proj, cellhd->zone);
268 else
269 sprintf(error_msg + strlen(error_msg),
270 "Location proj = %d (unknown), zone = %d\n",
271 loc_wind.proj, cellhd->zone);
272 }
273
274 if (proj_info != NULL) {
275 strcat(error_msg, _("Dataset PROJ_INFO is:\n"));
276 for (i_value = 0; i_value < proj_info->nitems; i_value++)
277 sprintf(error_msg + strlen(error_msg), "%s: %s\n",
278 proj_info->key[i_value],
279 proj_info->value[i_value]);
280 }
281 else {
282 strcat(error_msg, _("Dataset PROJ_INFO is:\n"));
283 if (cellhd->proj == PROJECTION_XY)
284 sprintf(error_msg + strlen(error_msg),
285 "Dataset proj = %d (unreferenced/unknown)\n",
286 cellhd->proj);
287 else if (cellhd->proj == PROJECTION_LL)
288 sprintf(error_msg + strlen(error_msg),
289 "Dataset proj = %d (lat/long)\n",
290 cellhd->proj);
291 else if (cellhd->proj == PROJECTION_UTM)
292 sprintf(error_msg + strlen(error_msg),
293 "Dataset proj = %d (UTM), zone = %d\n",
294 cellhd->proj, cellhd->zone);
295 else
296 sprintf(error_msg + strlen(error_msg),
297 "Dataset proj = %d (unknown), zone = %d\n",
298 cellhd->proj, cellhd->zone);
299 }
300 if (loc_wind.proj != cellhd->proj) {
301 strcat(error_msg, "\nDifference in: proj\n");
302 }
303 else {
304 strcat(error_msg, "\nDifference in: ");
305 switch (err) {
306 case -1:
307 strcat(error_msg, "proj\n");
308 break;
309 case -2:
310 strcat(error_msg, "units\n");
311 break;
312 case -3:
313 strcat(error_msg, "datum\n");
314 break;
315 case -4:
316 strcat(error_msg, "ellps, a, es\n");
317 break;
318 case -5:
319 strcat(error_msg, "zone\n");
320 break;
321 case -6:
322 strcat(error_msg, "south\n");
323 break;
324 case -7:
325 strcat(error_msg, "x_0\n");
326 break;
327 case -8:
328 strcat(error_msg, "y_0\n");
329 break;
330 case -9:
331 strcat(error_msg, "lon_0\n");
332 break;
333 case -10:
334 strcat(error_msg, "lat_0\n");
335 break;
336 case -11:
337 strcat(error_msg, "lat_1, lat2\n");
338 break;
339 }
340 }
341 }
342 else {
343 /* error in proj_units */
344 if (loc_proj_units != NULL) {
345 strcat(error_msg, "Location PROJ_UNITS is:\n");
346 for (i_value = 0; i_value < loc_proj_units->nitems;
347 i_value++)
348 sprintf(error_msg + strlen(error_msg), "%s: %s\n",
349 loc_proj_units->key[i_value],
350 loc_proj_units->value[i_value]);
351 strcat(error_msg, "\n");
352 }
353
354 if (proj_units != NULL) {
355 strcat(error_msg, "Dataset PROJ_UNITS is:\n");
356 for (i_value = 0; i_value < proj_units->nitems; i_value++)
357 sprintf(error_msg + strlen(error_msg), "%s: %s\n",
358 proj_units->key[i_value],
359 proj_units->value[i_value]);
360 }
361 }
362 if (!check_only) {
363 strcat(error_msg,
364 _("\nIn case of no significant differences in the projection definitions,"
365 " use the -o flag to ignore them and use"
366 " current location definition.\n"));
367 strcat(error_msg,
368 _("Consider generating a new location from the input dataset using "
369 "the 'location' parameter.\n"));
370 }
371
372 if (check_only)
373 msg_fn = G_message;
374 else
375 msg_fn = G_fatal_error;
376 msg_fn("%s", error_msg);
377 if (check_only) {
378 ds_close(hDS);
379 exit(EXIT_FAILURE);
380 }
381 }
382 else {
383 if (check_only)
384 msg_fn = G_message;
385 else
386 msg_fn = G_verbose_message;
387 msg_fn(_("Projection of input dataset and current location "
388 "appear to match"));
389 if (check_only) {
390 ds_close(hDS);
391 exit(EXIT_SUCCESS);
392 }
393 }
394 }
395 }
396