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