Line data    Source code

       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             : extern const MapInfoDatumInfo asDatumInfoList[];
      48             : extern const MapInfoSpheroidInfo asSpheroidInfoList[];
      49             : 
      50             : /************************************************************************/
      51             : /*                      MITABCoordSys2SpatialRef()                      */
      52             : /*                                                                      */
      53             : /*      Convert a MIF COORDSYS string into a new OGRSpatialReference    */
      54             : /*      object.                                                         */
      55             : /************************************************************************/
      56             : 
      57        7017 : OGRSpatialReference *MITABCoordSys2SpatialRef(const char *pszCoordSys)
      58             : 
      59             : {
      60             :     TABProjInfo sTABProj;
      61        7017 :     if (MITABCoordSys2TABProjInfo(pszCoordSys, &sTABProj) < 0)
      62        6855 :         return nullptr;
      63         162 :     OGRSpatialReference *poSR = TABFile::GetSpatialRefFromTABProj(sTABProj);
      64             : 
      65             :     // Report on translation.
      66         162 :     char *pszWKT = nullptr;
      67             : 
      68         162 :     poSR->exportToWkt(&pszWKT);
      69         162 :     if (pszWKT != nullptr)
      70             :     {
      71         162 :         CPLDebug("MITAB", "This CoordSys value:\n%s\nwas translated to:\n%s",
      72             :                  pszCoordSys, pszWKT);
      73         162 :         CPLFree(pszWKT);
      74             :     }
      75             : 
      76         162 :     return poSR;
      77             : }
      78             : 
      79             : /************************************************************************/
      80             : /*                      MITABSpatialRef2CoordSys()                      */
      81             : /*                                                                      */
      82             : /*      Converts a OGRSpatialReference object into a MIF COORDSYS       */
      83             : /*      string.                                                         */
      84             : /*                                                                      */
      85             : /*      The function returns a newly allocated string that should be    */
      86             : /*      CPLFree()'d by the caller.                                      */
      87             : /************************************************************************/
      88             : 
      89          81 : char *MITABSpatialRef2CoordSys(const OGRSpatialReference *poSR)
      90             : 
      91             : {
      92          81 :     if (poSR == nullptr)
      93           0 :         return nullptr;
      94             : 
      95             :     TABProjInfo sTABProj;
      96          81 :     int nParamCount = 0;
      97          81 :     TABFile::GetTABProjFromSpatialRef(poSR, sTABProj, nParamCount);
      98             : 
      99             :     // Do coordsys lookup.
     100          81 :     double dXMin = 0.0;
     101          81 :     double dYMin = 0.0;
     102          81 :     double dXMax = 0.0;
     103          81 :     double dYMax = 0.0;
     104          81 :     bool bHasBounds = false;
     105         136 :     if (sTABProj.nProjId > 1 &&
     106          55 :         MITABLookupCoordSysBounds(&sTABProj, dXMin, dYMin, dXMax, dYMax, true))
     107             :     {
     108           5 :         bHasBounds = true;
     109             :     }
     110             : 
     111             :     // Translate the units.
     112          81 :     const char *pszMIFUnits = TABUnitIdToString(sTABProj.nUnitsId);
     113             : 
     114             :     // Build coordinate system definition.
     115         162 :     CPLString osCoordSys;
     116             : 
     117          81 :     if (sTABProj.nProjId != 0)
     118             :     {
     119          66 :         osCoordSys.Printf("Earth Projection %d", sTABProj.nProjId);
     120             :     }
     121             :     else
     122             :     {
     123          15 :         osCoordSys.Printf("NonEarth Units");
     124             :     }
     125             : 
     126             :     // Append Datum.
     127          81 :     if (sTABProj.nProjId != 0)
     128             :     {
     129          66 :         osCoordSys += CPLSPrintf(", %d", sTABProj.nDatumId);
     130             : 
     131          66 :         if (sTABProj.nDatumId == 999 || sTABProj.nDatumId == 9999)
     132             :         {
     133             :             osCoordSys +=
     134           4 :                 CPLSPrintf(", %d, %.15g, %.15g, %.15g", sTABProj.nEllipsoidId,
     135             :                            sTABProj.dDatumShiftX, sTABProj.dDatumShiftY,
     136           4 :                            sTABProj.dDatumShiftZ);
     137             :         }
     138             : 
     139          66 :         if (sTABProj.nDatumId == 9999)
     140             :         {
     141             :             osCoordSys +=
     142             :                 CPLSPrintf(", %.15g, %.15g, %.15g, %.15g, %.15g",
     143             :                            sTABProj.adDatumParams[0], sTABProj.adDatumParams[1],
     144             :                            sTABProj.adDatumParams[2], sTABProj.adDatumParams[3],
     145           3 :                            sTABProj.adDatumParams[4]);
     146             :         }
     147             :     }
     148             : 
     149             :     // Append units.
     150          81 :     if (sTABProj.nProjId != 1 && pszMIFUnits != nullptr)
     151             :     {
     152          70 :         if (sTABProj.nProjId != 0)
     153          55 :             osCoordSys += ",";
     154             : 
     155          70 :         osCoordSys += CPLSPrintf(" \"%s\"", pszMIFUnits);
     156             :     }
     157             : 
     158             :     // Append Projection Params.
     159         321 :     for (int iParam = 0; iParam < nParamCount; iParam++)
     160         240 :         osCoordSys += CPLSPrintf(", %.15g", sTABProj.adProjParams[iParam]);
     161             : 
     162             :     // Append user bounds.
     163          81 :     if (bHasBounds)
     164             :     {
     165           5 :         if (fabs(dXMin - floor(dXMin + 0.5)) < 1e-8 &&
     166           5 :             fabs(dYMin - floor(dYMin + 0.5)) < 1e-8 &&
     167           5 :             fabs(dXMax - floor(dXMax + 0.5)) < 1e-8 &&
     168           5 :             fabs(dYMax - floor(dYMax + 0.5)) < 1e-8)
     169             :         {
     170           5 :             osCoordSys +=
     171             :                 CPLSPrintf(" Bounds (%d, %d) (%d, %d)", static_cast<int>(dXMin),
     172             :                            static_cast<int>(dYMin), static_cast<int>(dXMax),
     173           5 :                            static_cast<int>(dYMax));
     174             :         }
     175             :         else
     176             :         {
     177             :             osCoordSys += CPLSPrintf(" Bounds (%f, %f) (%f, %f)", dXMin, dYMin,
     178           0 :                                      dXMax, dYMax);
     179             :         }
     180             :     }
     181             : 
     182             :     // Report on translation.
     183          81 :     char *pszWKT = nullptr;
     184             : 
     185          81 :     poSR->exportToWkt(&pszWKT);
     186          81 :     if (pszWKT != nullptr)
     187             :     {
     188          81 :         CPLDebug("MITAB", "This WKT Projection:\n%s\n\ntranslates to:\n%s",
     189             :                  pszWKT, osCoordSys.c_str());
     190          81 :         CPLFree(pszWKT);
     191             :     }
     192             : 
     193          81 :     return CPLStrdup(osCoordSys.c_str());
     194             : }
     195             : 
     196             : /************************************************************************/
     197             : /*                      MITABExtractCoordSysBounds                      */
     198             : /*                                                                      */
     199             : /* Return true if MIF coordsys string contains a BOUNDS parameter and   */
     200             : /* Set x/y min/max values.                                              */
     201             : /************************************************************************/
     202             : 
     203          12 : bool MITABExtractCoordSysBounds(const char *pszCoordSys, double &dXMin,
     204             :                                 double &dYMin, double &dXMax, double &dYMax)
     205             : 
     206             : {
     207          12 :     if (pszCoordSys == nullptr)
     208           0 :         return false;
     209             : 
     210             :     char **papszFields =
     211          12 :         CSLTokenizeStringComplex(pszCoordSys, " ,()", TRUE, FALSE);
     212             : 
     213          12 :     int iBounds = CSLFindString(papszFields, "Bounds");
     214             : 
     215          12 :     if (iBounds >= 0 && iBounds + 4 < CSLCount(papszFields))
     216             :     {
     217          12 :         dXMin = CPLAtof(papszFields[++iBounds]);
     218          12 :         dYMin = CPLAtof(papszFields[++iBounds]);
     219          12 :         dXMax = CPLAtof(papszFields[++iBounds]);
     220          12 :         dYMax = CPLAtof(papszFields[++iBounds]);
     221          12 :         CSLDestroy(papszFields);
     222          12 :         return true;
     223             :     }
     224             : 
     225           0 :     CSLDestroy(papszFields);
     226           0 :     return false;
     227             : }
     228             : 
     229             : /**********************************************************************
     230             :  *                     MITABCoordSys2TABProjInfo()
     231             :  *
     232             :  * Convert a MIF COORDSYS string into a TABProjInfo structure.
     233             :  *
     234             :  * Returns 0 on success, -1 on error.
     235             :  **********************************************************************/
     236        7035 : int MITABCoordSys2TABProjInfo(const char *pszCoordSys, TABProjInfo *psProj)
     237             : 
     238             : {
     239             :     // Set all fields to zero, equivalent of NonEarth Units "mi"
     240        7035 :     memset(psProj, 0, sizeof(TABProjInfo));
     241             : 
     242        7035 :     if (pszCoordSys == nullptr)
     243        6855 :         return -1;
     244             : 
     245             :     // Parse the passed string into words.
     246         214 :     while (*pszCoordSys == ' ')
     247          34 :         pszCoordSys++;  // Eat leading spaces.
     248         180 :     if (STARTS_WITH_CI(pszCoordSys, "CoordSys") && pszCoordSys[8] != '\0')
     249          33 :         pszCoordSys += 9;
     250             : 
     251             :     char **papszFields =
     252         180 :         CSLTokenizeStringComplex(pszCoordSys, " ,", TRUE, FALSE);
     253             : 
     254             :     // Clip off Bounds information.
     255         180 :     int iBounds = CSLFindString(papszFields, "Bounds");
     256             : 
     257         316 :     while (iBounds != -1 && papszFields[iBounds] != nullptr)
     258             :     {
     259         136 :         CPLFree(papszFields[iBounds]);
     260         136 :         papszFields[iBounds] = nullptr;
     261         136 :         iBounds++;
     262             :     }
     263             : 
     264             :     // Fetch the projection.
     265         180 :     char **papszNextField = nullptr;
     266             : 
     267         327 :     if (CSLCount(papszFields) >= 3 && EQUAL(papszFields[0], "Earth") &&
     268         147 :         EQUAL(papszFields[1], "Projection"))
     269             :     {
     270         147 :         int nProjId = atoi(papszFields[2]);
     271         147 :         if (nProjId >= 3000)
     272           0 :             nProjId -= 3000;
     273         147 :         else if (nProjId >= 2000)
     274           0 :             nProjId -= 2000;
     275         147 :         else if (nProjId >= 1000)
     276           0 :             nProjId -= 1000;
     277             : 
     278         147 :         psProj->nProjId = static_cast<GByte>(nProjId);
     279         147 :         papszNextField = papszFields + 3;
     280             :     }
     281          33 :     else if (CSLCount(papszFields) >= 2 && EQUAL(papszFields[0], "NonEarth"))
     282             :     {
     283             :         // NonEarth Units "..." Bounds (x, y) (x, y)
     284          33 :         psProj->nProjId = 0;
     285          33 :         papszNextField = papszFields + 2;
     286             : 
     287          33 :         if (papszNextField[0] != nullptr && EQUAL(papszNextField[0], "Units"))
     288           0 :             papszNextField++;
     289             :     }
     290             :     else
     291             :     {
     292             :         // Invalid projection string ???
     293           0 :         if (CSLCount(papszFields) > 0)
     294           0 :             CPLError(CE_Warning, CPLE_IllegalArg,
     295             :                      "Failed parsing CoordSys: '%s'", pszCoordSys);
     296           0 :         CSLDestroy(papszFields);
     297           0 :         return -1;
     298             :     }
     299             : 
     300             :     // Fetch the datum information.
     301         180 :     int nDatum = 0;
     302             : 
     303         180 :     if (psProj->nProjId != 0 && CSLCount(papszNextField) > 0)
     304             :     {
     305         147 :         nDatum = atoi(papszNextField[0]);
     306         147 :         papszNextField++;
     307             :     }
     308             : 
     309         180 :     if ((nDatum == 999 || nDatum == 9999) && CSLCount(papszNextField) >= 4)
     310             :     {
     311           5 :         psProj->nEllipsoidId = static_cast<GByte>(atoi(papszNextField[0]));
     312           5 :         psProj->dDatumShiftX = CPLAtof(papszNextField[1]);
     313           5 :         psProj->dDatumShiftY = CPLAtof(papszNextField[2]);
     314           5 :         psProj->dDatumShiftZ = CPLAtof(papszNextField[3]);
     315           5 :         papszNextField += 4;
     316             : 
     317           5 :         if (nDatum == 9999 && CSLCount(papszNextField) >= 5)
     318             :         {
     319           3 :             psProj->adDatumParams[0] = CPLAtof(papszNextField[0]);
     320           3 :             psProj->adDatumParams[1] = CPLAtof(papszNextField[1]);
     321           3 :             psProj->adDatumParams[2] = CPLAtof(papszNextField[2]);
     322           3 :             psProj->adDatumParams[3] = CPLAtof(papszNextField[3]);
     323           3 :             psProj->adDatumParams[4] = CPLAtof(papszNextField[4]);
     324           3 :             papszNextField += 5;
     325             :         }
     326             :     }
     327         175 :     else if (nDatum != 999 && nDatum != 9999)
     328             :     {
     329             :         // Find the datum, and collect its parameters if possible.
     330         175 :         const MapInfoDatumInfo *psDatumInfo = nullptr;
     331             : 
     332         175 :         int iDatum = 0;  // Used after for.
     333        5794 :         for (; asDatumInfoList[iDatum].nMapInfoDatumID != -1; iDatum++)
     334             :         {
     335        5794 :             if (asDatumInfoList[iDatum].nMapInfoDatumID == nDatum)
     336             :             {
     337         175 :                 psDatumInfo = asDatumInfoList + iDatum;
     338         175 :                 break;
     339             :             }
     340             :         }
     341             : 
     342         175 :         if (asDatumInfoList[iDatum].nMapInfoDatumID == -1)
     343             :         {
     344             :             // Use WGS84.
     345           0 :             psDatumInfo = asDatumInfoList + 0;
     346             :         }
     347             : 
     348         175 :         if (psDatumInfo != nullptr)
     349             :         {
     350         175 :             psProj->nEllipsoidId = static_cast<GByte>(psDatumInfo->nEllipsoid);
     351         175 :             psProj->nDatumId =
     352         175 :                 static_cast<GInt16>(psDatumInfo->nMapInfoDatumID);
     353         175 :             psProj->dDatumShiftX = psDatumInfo->dfShiftX;
     354         175 :             psProj->dDatumShiftY = psDatumInfo->dfShiftY;
     355         175 :             psProj->dDatumShiftZ = psDatumInfo->dfShiftZ;
     356         175 :             psProj->adDatumParams[0] = psDatumInfo->dfDatumParm0;
     357         175 :             psProj->adDatumParams[1] = psDatumInfo->dfDatumParm1;
     358         175 :             psProj->adDatumParams[2] = psDatumInfo->dfDatumParm2;
     359         175 :             psProj->adDatumParams[3] = psDatumInfo->dfDatumParm3;
     360         175 :             psProj->adDatumParams[4] = psDatumInfo->dfDatumParm4;
     361             :         }
     362             :     }
     363             : 
     364             :     // Fetch the units string.
     365         180 :     if (CSLCount(papszNextField) > 0)
     366             :     {
     367         166 :         if (isdigit(static_cast<unsigned char>(papszNextField[0][0])))
     368             :         {
     369           0 :             psProj->nUnitsId = static_cast<GByte>(atoi(papszNextField[0]));
     370             :         }
     371             :         else
     372             :         {
     373         166 :             psProj->nUnitsId =
     374         166 :                 static_cast<GByte>(TABUnitIdFromString(papszNextField[0]));
     375             :         }
     376         166 :         papszNextField++;
     377             :     }
     378             : 
     379             :     // Finally the projection parameters.
     380         802 :     for (int iParam = 0; iParam < 7 && CSLCount(papszNextField) > 0; iParam++)
     381             :     {
     382         622 :         psProj->adProjParams[iParam] = CPLAtof(papszNextField[0]);
     383         622 :         papszNextField++;
     384             :     }
     385             : 
     386         180 :     CSLDestroy(papszFields);
     387             : 
     388         180 :     return 0;
     389             : }