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