1 /**********************************************************************
2 *
3 * Name: mitab_coordsys.cpp
4 * Project: MapInfo TAB Read/Write library
5 * Language: C++
6 * Purpose: Implementation translation between MIF CoordSys format, and
7 * and OGRSpatialRef format.
8 * Author: Frank Warmerdam, warmerdam@pobox.com
9 *
10 **********************************************************************
11 * Copyright (c) 1999-2001, Frank Warmerdam
12 *
13 * Permission is hereby granted, free of charge, to any person obtaining a
14 * copy of this software and associated documentation files (the "Software"),
15 * to deal in the Software without restriction, including without limitation
16 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 * and/or sell copies of the Software, and to permit persons to whom the
18 * Software is furnished to do so, subject to the following conditions:
19 *
20 * The above copyright notice and this permission notice shall be included
21 * in all copies or substantial portions of the Software.
22 *
23 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
24 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 * DEALINGS IN THE SOFTWARE.
30 **********************************************************************/
31
32 #include "cpl_port.h"
33 #include "mitab.h"
34 #include "mitab_utils.h"
35
36 #include <cmath>
37 #include <cstdlib>
38 #include <cstring>
39 #include <string>
40
41 #include "cpl_conv.h"
42 #include "cpl_error.h"
43 #include "cpl_string.h"
44 #include "mitab_priv.h"
45 #include "ogr_spatialref.h"
46
47 CPL_CVSID("$Id: mitab_coordsys.cpp fa752ad6eabafaf630a704e1892a9d837d683cb3 2021-03-06 17:04:38 +0100 Even Rouault $")
48
49 extern const MapInfoDatumInfo asDatumInfoList[];
50 extern const MapInfoSpheroidInfo asSpheroidInfoList[];
51
52 /************************************************************************/
53 /* MITABCoordSys2SpatialRef() */
54 /* */
55 /* Convert a MIF COORDSYS string into a new OGRSpatialReference */
56 /* object. */
57 /************************************************************************/
58
MITABCoordSys2SpatialRef(const char * pszCoordSys)59 OGRSpatialReference *MITABCoordSys2SpatialRef( const char * pszCoordSys )
60
61 {
62 TABProjInfo sTABProj;
63 if(MITABCoordSys2TABProjInfo(pszCoordSys, &sTABProj) < 0 )
64 return nullptr;
65 OGRSpatialReference *poSR = TABFile::GetSpatialRefFromTABProj(sTABProj);
66
67 // Report on translation.
68 char *pszWKT = nullptr;
69
70 poSR->exportToWkt(&pszWKT);
71 if( pszWKT != nullptr )
72 {
73 CPLDebug("MITAB",
74 "This CoordSys value:\n%s\nwas translated to:\n%s",
75 pszCoordSys, pszWKT);
76 CPLFree(pszWKT);
77 }
78
79 return poSR;
80 }
81
82 /************************************************************************/
83 /* MITABSpatialRef2CoordSys() */
84 /* */
85 /* Converts a OGRSpatialReference object into a MIF COORDSYS */
86 /* string. */
87 /* */
88 /* The function returns a newly allocated string that should be */
89 /* CPLFree()'d by the caller. */
90 /************************************************************************/
91
MITABSpatialRef2CoordSys(const OGRSpatialReference * poSR)92 char *MITABSpatialRef2CoordSys( const OGRSpatialReference * poSR )
93
94 {
95 if( poSR == nullptr )
96 return nullptr;
97
98 TABProjInfo sTABProj;
99 int nParamCount = 0;
100 TABFile::GetTABProjFromSpatialRef(poSR, sTABProj, nParamCount);
101
102 // Do coordsys lookup.
103 double dXMin = 0.0;
104 double dYMin = 0.0;
105 double dXMax = 0.0;
106 double dYMax = 0.0;
107 bool bHasBounds = false;
108 if( sTABProj.nProjId > 1 &&
109 MITABLookupCoordSysBounds(&sTABProj,
110 dXMin, dYMin,
111 dXMax, dYMax, true) )
112 {
113 bHasBounds = true;
114 }
115
116 // Translate the units.
117 const char *pszMIFUnits = TABUnitIdToString(sTABProj.nUnitsId);
118
119 // Build coordinate system definition.
120 CPLString osCoordSys;
121
122 if( sTABProj.nProjId != 0 )
123 {
124 osCoordSys.Printf("Earth Projection %d", sTABProj.nProjId);
125 }
126 else
127 {
128 osCoordSys.Printf("NonEarth Units");
129 }
130
131 // Append Datum.
132 if( sTABProj.nProjId != 0 )
133 {
134 osCoordSys += CPLSPrintf(", %d", sTABProj.nDatumId);
135
136 if( sTABProj.nDatumId == 999 || sTABProj.nDatumId == 9999 )
137 {
138 osCoordSys +=
139 CPLSPrintf(", %d, %.15g, %.15g, %.15g", sTABProj.nEllipsoidId,
140 sTABProj.dDatumShiftX, sTABProj.dDatumShiftY,
141 sTABProj.dDatumShiftZ);
142 }
143
144 if( sTABProj.nDatumId == 9999 )
145 {
146 osCoordSys +=
147 CPLSPrintf(", %.15g, %.15g, %.15g, %.15g, %.15g",
148 sTABProj.adDatumParams[0], sTABProj.adDatumParams[1],
149 sTABProj.adDatumParams[2], sTABProj.adDatumParams[3],
150 sTABProj.adDatumParams[4]);
151 }
152 }
153
154 // Append units.
155 if( sTABProj.nProjId != 1 && pszMIFUnits != nullptr )
156 {
157 if( sTABProj.nProjId != 0 )
158 osCoordSys += ",";
159
160 osCoordSys += CPLSPrintf(" \"%s\"", pszMIFUnits);
161 }
162
163 // Append Projection Params.
164 for( int iParam = 0; iParam < nParamCount; iParam++ )
165 osCoordSys += CPLSPrintf(", %.15g", sTABProj.adProjParams[iParam]);
166
167 // Append user bounds.
168 if( bHasBounds )
169 {
170 if( fabs(dXMin - floor(dXMin + 0.5)) < 1e-8 &&
171 fabs(dYMin - floor(dYMin + 0.5)) < 1e-8 &&
172 fabs(dXMax - floor(dXMax + 0.5)) < 1e-8 &&
173 fabs(dYMax - floor(dYMax + 0.5)) < 1e-8 )
174 {
175 osCoordSys +=
176 CPLSPrintf(" Bounds (%d, %d) (%d, %d)",
177 static_cast<int>(dXMin), static_cast<int>(dYMin),
178 static_cast<int>(dXMax), static_cast<int>(dYMax));
179 }
180 else
181 {
182 osCoordSys += CPLSPrintf(" Bounds (%f, %f) (%f, %f)",
183 dXMin, dYMin, dXMax, dYMax);
184 }
185 }
186
187 // Report on translation.
188 char *pszWKT = nullptr;
189
190 poSR->exportToWkt(&pszWKT);
191 if( pszWKT != nullptr )
192 {
193 CPLDebug("MITAB",
194 "This WKT Projection:\n%s\n\ntranslates to:\n%s",
195 pszWKT, osCoordSys.c_str());
196 CPLFree(pszWKT);
197 }
198
199 return CPLStrdup(osCoordSys.c_str());
200 }
201
202 /************************************************************************/
203 /* MITABExtractCoordSysBounds */
204 /* */
205 /* Return true if MIF coordsys string contains a BOUNDS parameter and */
206 /* Set x/y min/max values. */
207 /************************************************************************/
208
MITABExtractCoordSysBounds(const char * pszCoordSys,double & dXMin,double & dYMin,double & dXMax,double & dYMax)209 bool MITABExtractCoordSysBounds( const char * pszCoordSys,
210 double &dXMin, double &dYMin,
211 double &dXMax, double &dYMax )
212
213 {
214 if( pszCoordSys == nullptr )
215 return false;
216
217 char **papszFields =
218 CSLTokenizeStringComplex(pszCoordSys, " ,()", TRUE, FALSE);
219
220 int iBounds = CSLFindString(papszFields, "Bounds");
221
222 if (iBounds >= 0 && iBounds + 4 < CSLCount(papszFields))
223 {
224 dXMin = CPLAtof(papszFields[++iBounds]);
225 dYMin = CPLAtof(papszFields[++iBounds]);
226 dXMax = CPLAtof(papszFields[++iBounds]);
227 dYMax = CPLAtof(papszFields[++iBounds]);
228 CSLDestroy(papszFields);
229 return true;
230 }
231
232 CSLDestroy(papszFields);
233 return false;
234 }
235
236 /**********************************************************************
237 * MITABCoordSys2TABProjInfo()
238 *
239 * Convert a MIF COORDSYS string into a TABProjInfo structure.
240 *
241 * Returns 0 on success, -1 on error.
242 **********************************************************************/
MITABCoordSys2TABProjInfo(const char * pszCoordSys,TABProjInfo * psProj)243 int MITABCoordSys2TABProjInfo(const char * pszCoordSys, TABProjInfo *psProj)
244
245 {
246 // Set all fields to zero, equivalent of NonEarth Units "mi"
247 memset(psProj, 0, sizeof(TABProjInfo));
248
249 if( pszCoordSys == nullptr )
250 return -1;
251
252 // Parse the passed string into words.
253 while(*pszCoordSys == ' ')
254 pszCoordSys++; // Eat leading spaces.
255 if( STARTS_WITH_CI(pszCoordSys, "CoordSys") && pszCoordSys[8] != '\0' )
256 pszCoordSys += 9;
257
258 char **papszFields =
259 CSLTokenizeStringComplex(pszCoordSys, " ,", TRUE, FALSE);
260
261 // Clip off Bounds information.
262 int iBounds = CSLFindString(papszFields, "Bounds");
263
264 while( iBounds != -1 && papszFields[iBounds] != nullptr )
265 {
266 CPLFree( papszFields[iBounds] );
267 papszFields[iBounds] = nullptr;
268 iBounds++;
269 }
270
271 // Fetch the projection.
272 char **papszNextField = nullptr;
273
274 if( CSLCount(papszFields) >= 3 &&
275 EQUAL(papszFields[0], "Earth") &&
276 EQUAL(papszFields[1], "Projection") )
277 {
278 int nProjId = atoi(papszFields[2]);
279 if (nProjId >= 3000) nProjId -= 3000;
280 else if (nProjId >= 2000) nProjId -= 2000;
281 else if (nProjId >= 1000) nProjId -= 1000;
282
283 psProj->nProjId = static_cast<GByte>(nProjId);
284 papszNextField = papszFields + 3;
285 }
286 else if (CSLCount(papszFields) >= 2 &&
287 EQUAL(papszFields[0],"NonEarth") )
288 {
289 // NonEarth Units "..." Bounds (x, y) (x, y)
290 psProj->nProjId = 0;
291 papszNextField = papszFields + 2;
292
293 if( papszNextField[0] != nullptr && EQUAL(papszNextField[0], "Units") )
294 papszNextField++;
295 }
296 else
297 {
298 // Invalid projection string ???
299 if (CSLCount(papszFields) > 0)
300 CPLError(CE_Warning, CPLE_IllegalArg,
301 "Failed parsing CoordSys: '%s'", pszCoordSys);
302 CSLDestroy(papszFields);
303 return -1;
304 }
305
306 // Fetch the datum information.
307 int nDatum = 0;
308
309 if( psProj->nProjId != 0 && CSLCount(papszNextField) > 0 )
310 {
311 nDatum = atoi(papszNextField[0]);
312 papszNextField++;
313 }
314
315 if( (nDatum == 999 || nDatum == 9999) &&
316 CSLCount(papszNextField) >= 4 )
317 {
318 psProj->nEllipsoidId = static_cast<GByte>(atoi(papszNextField[0]));
319 psProj->dDatumShiftX = CPLAtof(papszNextField[1]);
320 psProj->dDatumShiftY = CPLAtof(papszNextField[2]);
321 psProj->dDatumShiftZ = CPLAtof(papszNextField[3]);
322 papszNextField += 4;
323
324 if( nDatum == 9999 &&
325 CSLCount(papszNextField) >= 5 )
326 {
327 psProj->adDatumParams[0] = CPLAtof(papszNextField[0]);
328 psProj->adDatumParams[1] = CPLAtof(papszNextField[1]);
329 psProj->adDatumParams[2] = CPLAtof(papszNextField[2]);
330 psProj->adDatumParams[3] = CPLAtof(papszNextField[3]);
331 psProj->adDatumParams[4] = CPLAtof(papszNextField[4]);
332 papszNextField += 5;
333 }
334 }
335 else if (nDatum != 999 && nDatum != 9999)
336 {
337 // Find the datum, and collect its parameters if possible.
338 const MapInfoDatumInfo *psDatumInfo = nullptr;
339
340 int iDatum = 0; // Used after for.
341 for( ; asDatumInfoList[iDatum].nMapInfoDatumID != -1; iDatum++ )
342 {
343 if( asDatumInfoList[iDatum].nMapInfoDatumID == nDatum )
344 {
345 psDatumInfo = asDatumInfoList + iDatum;
346 break;
347 }
348 }
349
350 if( asDatumInfoList[iDatum].nMapInfoDatumID == -1 )
351 {
352 // Use WGS84.
353 psDatumInfo = asDatumInfoList + 0;
354 }
355
356 if( psDatumInfo != nullptr )
357 {
358 psProj->nEllipsoidId = static_cast<GByte>(psDatumInfo->nEllipsoid);
359 psProj->nDatumId =
360 static_cast<GInt16>(psDatumInfo->nMapInfoDatumID);
361 psProj->dDatumShiftX = psDatumInfo->dfShiftX;
362 psProj->dDatumShiftY = psDatumInfo->dfShiftY;
363 psProj->dDatumShiftZ = psDatumInfo->dfShiftZ;
364 psProj->adDatumParams[0] = psDatumInfo->dfDatumParm0;
365 psProj->adDatumParams[1] = psDatumInfo->dfDatumParm1;
366 psProj->adDatumParams[2] = psDatumInfo->dfDatumParm2;
367 psProj->adDatumParams[3] = psDatumInfo->dfDatumParm3;
368 psProj->adDatumParams[4] = psDatumInfo->dfDatumParm4;
369 }
370 }
371
372 // Fetch the units string.
373 if( CSLCount(papszNextField) > 0 )
374 {
375 if( isdigit(papszNextField[0][0]) )
376 {
377 psProj->nUnitsId =
378 static_cast<GByte>(atoi(papszNextField[0]));
379 }
380 else
381 {
382 psProj->nUnitsId =
383 static_cast<GByte>(TABUnitIdFromString(papszNextField[0]));
384 }
385 papszNextField++;
386 }
387
388 // Finally the projection parameters.
389 for(int iParam = 0; iParam < 6 && CSLCount(papszNextField) > 0; iParam++)
390 {
391 psProj->adProjParams[iParam] = CPLAtof(papszNextField[0]);
392 papszNextField++;
393 }
394
395 CSLDestroy(papszFields);
396
397 return 0;
398 }
399