Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  S-101 driver
       4             :  * Purpose:  Implements OGRS101Reader
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2026, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "ogr_s101.h"
      14             : #include "ogrs101readerconstants.h"
      15             : 
      16             : #include <memory>
      17             : 
      18             : /************************************************************************/
      19             : /*                           LaunderCRSName()                           */
      20             : /************************************************************************/
      21             : 
      22             : /** Launder a CRS name to be layer name friendly. */
      23             : /* static */ std::string
      24          57 : OGRS101Reader::LaunderCRSName(const OGRSpatialReference &oSRS)
      25             : {
      26          57 :     const char *pszSRSName = oSRS.GetName();
      27         114 :     return CPLString(pszSRSName ? pszSRSName : "(null)")
      28         114 :         .replaceAll("WGS 84 + ", "")
      29         114 :         .replaceAll(" depth", "");
      30             : }
      31             : 
      32             : /************************************************************************/
      33             : /*                              ReadCSID()                              */
      34             : /************************************************************************/
      35             : 
      36             : /** Read the (mandatory) CRS record.
      37             :  *
      38             :  * The horizontal CRS is mandatory and WGS 84 geographic.
      39             :  * There might be 0..n vertical CRS.
      40             :  */
      41         451 : bool OGRS101Reader::ReadCSID(const DDFRecord *poRecord)
      42             : {
      43             : 
      44         451 :     const auto poField = poRecord->FindField(CSID_FIELD);
      45         451 :     if (!poField)
      46           1 :         return EMIT_ERROR("CSID field not found");
      47             : 
      48             :     // Record name
      49             :     const RecordName nRCNM =
      50         450 :         poRecord->GetIntSubfield(poField, RCNM_SUBFIELD, 0);
      51         452 :     if (nRCNM != RECORD_NAME_CRS &&
      52           2 :         !EMIT_ERROR_OR_WARNING("Invalid value for RCNM subfield of CSID."))
      53           1 :         return false;
      54             : 
      55             :     // Record identifier
      56         449 :     const int nRCID = poRecord->GetIntSubfield(poField, RCID_SUBFIELD, 0);
      57             :     // Only one CRS record expected
      58         451 :     if (nRCID != 1 &&
      59           2 :         !EMIT_ERROR_OR_WARNING("Invalid value for RCID subfield of CSID."))
      60           1 :         return false;
      61             : 
      62             :     // Number of CRS Components
      63         448 :     constexpr const char *NCRC_SUBFIELD = "NCRC";
      64         448 :     int nNCRC = poRecord->GetIntSubfield(poField, NCRC_SUBFIELD, 0);
      65         450 :     if (nNCRC < 1 &&
      66           2 :         !EMIT_ERROR_OR_WARNING("Invalid value for NCRC subfield of CSID."))
      67           1 :         return false;
      68             : 
      69             :     // Fields in that record have normally the following sequence:
      70             :     // - CSID
      71             :     // - CRSH: horizontal CRS
      72             :     // - CRSH: vertical CRS 1
      73             :     //      - CSAX: coordinate system axis
      74             :     //      - VDAT: vertical datum
      75             :     // ...
      76             :     // - CRSH: vertical CRS N
      77             :     //      - CSAX: coordinate system axis
      78             :     //      - VDAT: vertical datum
      79             :     // Establish a mapping from the CRSH instance count to the CSAX/VDAT
      80             :     // instance count
      81             : 
      82         447 :     const int nFieldCount = poRecord->GetFieldCount();
      83         894 :     std::vector<int> anCSAXInstanceForCRSH;
      84         894 :     std::vector<int> anVDATInstanceForCRSH;
      85         447 :     int iInstanceCRSH = -1;
      86         447 :     int iInstanceCSAX = -1;
      87         447 :     int iInstanceVDAT = -1;
      88        1182 :     for (int i = 1; i < nFieldCount; ++i)
      89             :     {
      90             :         const char *pszFieldName =
      91         739 :             poRecord->GetField(i)->GetFieldDefn()->GetName();
      92         739 :         if (strcmp(pszFieldName, CRSH_FIELD) == 0)
      93             :         {
      94         547 :             ++iInstanceCRSH;
      95         547 :             anCSAXInstanceForCRSH.push_back(-1);
      96         547 :             anVDATInstanceForCRSH.push_back(-1);
      97             :         }
      98         192 :         else if (strcmp(pszFieldName, CSAX_FIELD) == 0)
      99             :         {
     100         100 :             ++iInstanceCSAX;
     101         100 :             if (iInstanceCRSH >= 0)
     102             :             {
     103          98 :                 CPLAssert(static_cast<size_t>(iInstanceCRSH) <
     104             :                           anCSAXInstanceForCRSH.size());
     105          98 :                 if (anCSAXInstanceForCRSH[iInstanceCRSH] >= 0)
     106             :                 {
     107           4 :                     if (iInstanceCRSH != 0 &&
     108           2 :                         !EMIT_ERROR_OR_WARNING(
     109             :                             CPLSPrintf("Several CSAX fields associated to CRSH "
     110             :                                        "instance of index %d.",
     111             :                                        iInstanceCRSH)))
     112             :                     {
     113           1 :                         return false;
     114             :                     }
     115             :                 }
     116             :                 else
     117             :                 {
     118          96 :                     anCSAXInstanceForCRSH[iInstanceCRSH] = iInstanceCSAX;
     119             :                 }
     120             :             }
     121             :             else
     122             :             {
     123           2 :                 if (!EMIT_ERROR_OR_WARNING("CSAX field found before CRSH."))
     124             :                 {
     125           1 :                     return false;
     126             :                 }
     127             :             }
     128             :         }
     129          92 :         else if (strcmp(pszFieldName, VDAT_FIELD) == 0)
     130             :         {
     131          92 :             ++iInstanceVDAT;
     132          92 :             if (iInstanceCRSH >= 0)
     133             :             {
     134          90 :                 CPLAssert(static_cast<size_t>(iInstanceCRSH) <
     135             :                           anVDATInstanceForCRSH.size());
     136          90 :                 if (anVDATInstanceForCRSH[iInstanceCRSH] >= 0)
     137             :                 {
     138           2 :                     if (!EMIT_ERROR_OR_WARNING(
     139             :                             CPLSPrintf("Several VDAT fields associated to CRSH "
     140             :                                        "instance of index %d.",
     141             :                                        iInstanceCRSH)))
     142             :                     {
     143           1 :                         return false;
     144             :                     }
     145             :                 }
     146             :                 else
     147             :                 {
     148          88 :                     anVDATInstanceForCRSH[iInstanceCRSH] = iInstanceVDAT;
     149             :                 }
     150             :             }
     151             :             else
     152             :             {
     153           2 :                 if (!EMIT_ERROR_OR_WARNING("VDAT field found before CRSH."))
     154             :                 {
     155           1 :                     return false;
     156             :                 }
     157             :             }
     158             :         }
     159           0 :         else if (!EMIT_ERROR_OR_WARNING(
     160             :                      CPLSPrintf("Unexpected field found in CRS record: %s.",
     161             :                                 pszFieldName)))
     162             :         {
     163           0 :             return false;
     164             :         }
     165             :     }
     166         448 :     if (nNCRC != iInstanceCRSH + 1 &&
     167           5 :         !EMIT_ERROR_OR_WARNING(
     168             :             "NCRC field of CSID is not consistent with number of CRSH fields."))
     169             :     {
     170           2 :         return false;
     171             :     }
     172             : 
     173         441 :     if (iInstanceCRSH < 0)
     174             :     {
     175           1 :         return EMIT_ERROR_OR_WARNING("No CRSH field.");
     176             :     }
     177             : 
     178         440 :     nNCRC = std::clamp(nNCRC, 1, iInstanceCRSH + 1);
     179             : 
     180             :     // Subfields of field CRSH
     181         440 :     constexpr const char *CRIX_SUBFIELD = "CRIX";
     182         440 :     constexpr const char *CRST_SUBFIELD = "CRST";
     183         440 :     constexpr const char *CSTY_SUBFIELD = "CSTY";
     184         440 :     constexpr const char *CRNM_SUBFIELD = "CRNM";
     185         440 :     constexpr const char *CRSI_SUBFIELD = "CRSI";
     186         440 :     constexpr const char *CRSS_SUBFIELD = "CRSS";
     187             : 
     188         880 :     const auto apoCRSHFields = poRecord->GetFields(CRSH_FIELD);
     189             : 
     190             :     // Read horizontal CRS definition (required to be WGS 84 by S-101 spec)
     191             :     {
     192         440 :         CPLAssert(!apoCRSHFields.empty());
     193         440 :         const auto poCRSHField = apoCRSHFields[0];
     194             : 
     195             :         // CRS Index
     196             :         const int nCRIX =
     197         440 :             poRecord->GetIntSubfield(poCRSHField, CRIX_SUBFIELD, 0);
     198             :         // 1 for the horizontal CRS
     199         440 :         if (nCRIX != 1 && !EMIT_ERROR_OR_WARNING(
     200             :                               "Invalid value for CRIX field of CRSH idx 0."))
     201           1 :             return false;
     202             : 
     203             :         // CRS Type
     204             :         const int nCRST =
     205         439 :             poRecord->GetIntSubfield(poCRSHField, CRST_SUBFIELD, 0);
     206         439 :         constexpr int CRS_TYPE_GEOGRAPHIC_2D = 1;
     207         441 :         if (nCRST != CRS_TYPE_GEOGRAPHIC_2D &&
     208           2 :             !EMIT_ERROR_OR_WARNING(
     209             :                 "Invalid value for CRST field of CRSH idx 0."))
     210           1 :             return false;
     211             : 
     212             :         // Coordinate System Type
     213             :         const int nCSTY =
     214         438 :             poRecord->GetIntSubfield(poCRSHField, CSTY_SUBFIELD, 0);
     215         438 :         constexpr int CS_TYPE_ELLIPSOIDAL = 1;
     216         440 :         if (nCSTY != CS_TYPE_ELLIPSOIDAL &&
     217           2 :             !EMIT_ERROR_OR_WARNING(
     218             :                 "Invalid value for CSTY field of CRSH idx 0."))
     219           1 :             return false;
     220             : 
     221             :         // CRS Name
     222             :         const char *pszCRNM =
     223         437 :             poRecord->GetStringSubfield(poCRSHField, CRNM_SUBFIELD, 0);
     224         437 :         if (!pszCRNM || strcmp(pszCRNM, "WGS84") != 0)
     225             :         {
     226           2 :             if (!EMIT_ERROR_OR_WARNING(
     227             :                     "Invalid value for CRNM field of CRSH idx 0."))
     228           1 :                 return false;
     229             :         }
     230             : 
     231             :         // CRS Identifier
     232             :         const char *pszCRSI =
     233         436 :             poRecord->GetStringSubfield(poCRSHField, CRSI_SUBFIELD, 0);
     234         436 :         if (!pszCRSI || strcmp(pszCRSI, "4326") != 0)
     235             :         {
     236           2 :             if (!EMIT_ERROR_OR_WARNING(
     237             :                     "Invalid value for CRSI field of CRSH idx 0."))
     238           1 :                 return false;
     239             :         }
     240             : 
     241             :         // CRS Source
     242             :         const int nCRSS =
     243         435 :             poRecord->GetIntSubfield(poCRSHField, CRSS_SUBFIELD, 0);
     244         435 :         constexpr int CRS_SOURCE_EPSG = 2;
     245         437 :         if (nCRSS != CRS_SOURCE_EPSG &&
     246           2 :             !EMIT_ERROR_OR_WARNING(
     247             :                 "Invalid value for CRSS field of CRSH idx 0."))
     248           1 :             return false;
     249             : 
     250         434 :         CPLAssert(!anCSAXInstanceForCRSH.empty());
     251         434 :         const int nCSAXIdx = anCSAXInstanceForCRSH[0];
     252         436 :         if (nCSAXIdx >= 0 &&
     253           2 :             !EMIT_ERROR_OR_WARNING(
     254             :                 "Unexpected CSAX field associated to first CRSH field."))
     255           1 :             return false;
     256             : 
     257         433 :         CPLAssert(!anVDATInstanceForCRSH.empty());
     258         433 :         const int nVDATIdx = anVDATInstanceForCRSH[0];
     259         435 :         if (nVDATIdx >= 0 &&
     260           2 :             !EMIT_ERROR_OR_WARNING(
     261             :                 "Unexpected VDAT field associated to first CRSH field."))
     262           1 :             return false;
     263             : 
     264         432 :         const int nEPSGCode = pszCRSI ? atoi(pszCRSI) : 0;
     265         432 :         if (nEPSGCode > 0)
     266             :         {
     267         432 :             OGRSpatialReference oSRSHorizontal;
     268         432 :             oSRSHorizontal.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     269         433 :             if (oSRSHorizontal.importFromEPSG(nEPSGCode) != OGRERR_NONE &&
     270           1 :                 m_bStrict)
     271           0 :                 return false;
     272             : 
     273         432 :             m_oMapSRS[nCRIX] = std::move(oSRSHorizontal);
     274             :         }
     275             :         else
     276             :         {
     277           0 :             return EMIT_ERROR(
     278             :                 CPLSPrintf("Invalid CRS EPSG code: %d.", nEPSGCode));
     279             :         }
     280             :     }
     281             : 
     282         432 :     const auto &oSRSHorizontal = m_oMapSRS[HORIZONTAL_CRS_ID];
     283             : 
     284             :     // Loop through optional vertical CRS definitions
     285         523 :     for (int iCRSHIdx = 1; iCRSHIdx < nNCRC; ++iCRSHIdx)
     286             :     {
     287         102 :         const auto poCRSHField = apoCRSHFields[iCRSHIdx];
     288             : 
     289             :         // CRS Index
     290             :         const int nCRIX =
     291         102 :             poRecord->GetIntSubfield(poCRSHField, CRIX_SUBFIELD, 0);
     292         104 :         if (!(nCRIX >= 2 && nCRIX <= nNCRC) &&
     293           2 :             !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     294             :                 "Invalid value for CRIX field of CRSH idx %d.", iCRSHIdx)))
     295             :         {
     296          11 :             return false;
     297             :         }
     298             : 
     299         101 :         if (cpl::contains(m_oMapSRS, nCRIX) &&
     300           0 :             !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     301             :                 "Several CRSH field instances with same CRIX = %d.", iCRSHIdx)))
     302             :         {
     303           0 :             return false;
     304             :         }
     305             : 
     306             :         // CRS Type
     307             :         const int nCRST =
     308         101 :             poRecord->GetIntSubfield(poCRSHField, CRST_SUBFIELD, 0);
     309         101 :         constexpr int CRS_TYPE_VERTICAL = 5;
     310         103 :         if (nCRST != CRS_TYPE_VERTICAL &&
     311           2 :             !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     312             :                 "Invalid value for CRST field of CRSH idx %d.", iCRSHIdx)))
     313             :         {
     314           1 :             return false;
     315             :         }
     316             : 
     317             :         // Coordinate System Type
     318             :         const int nCSTY =
     319         100 :             poRecord->GetIntSubfield(poCRSHField, CSTY_SUBFIELD, 0);
     320         100 :         constexpr int CS_TYPE_VERTICAL = 3;
     321         102 :         if (nCSTY != CS_TYPE_VERTICAL &&
     322           2 :             !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     323             :                 "Invalid value for CSTY field of CRSH idx %d.", iCRSHIdx)))
     324           1 :             return false;
     325             : 
     326             :         // CRS Name
     327             :         const char *pszCRNM =
     328          99 :             poRecord->GetStringSubfield(poCRSHField, CRNM_SUBFIELD, 0);
     329          99 :         if (!pszCRNM &&
     330           0 :             !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     331             :                 "No value for CRNM field of CRSH idx %d.", iCRSHIdx)))
     332           0 :             return false;
     333          99 :         if (!pszCRNM)
     334           0 :             pszCRNM = "(null)";
     335             : 
     336             :         // CRS Source
     337             :         const int nCRSS =
     338          99 :             poRecord->GetIntSubfield(poCRSHField, CRSS_SUBFIELD, 0);
     339          99 :         constexpr int CRS_SOURCE_NOT_APPLICABLE = 255;
     340         101 :         if (nCRSS != CRS_SOURCE_NOT_APPLICABLE &&
     341           2 :             !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     342             :                 "Invalid value for CRSS field of CRSH idx %d.", iCRSHIdx)))
     343           1 :             return false;
     344             : 
     345             :         // Read associated CSAX field
     346          98 :         CPLAssert(static_cast<size_t>(iCRSHIdx) < anCSAXInstanceForCRSH.size());
     347          98 :         const int nCSAXIdx = anCSAXInstanceForCRSH[iCRSHIdx];
     348          98 :         if (nCSAXIdx < 0 && !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     349             :                                 "No CSAX field for CRSH idx %d.", iCRSHIdx)))
     350             :         {
     351           1 :             return false;
     352             :         }
     353             : 
     354          97 :         if (nCSAXIdx >= 0)
     355             :         {
     356          92 :             const auto poCSAXField = poRecord->FindField(CSAX_FIELD, nCSAXIdx);
     357             : 
     358             :             // Axis type
     359          92 :             const int nAXTY = poRecord->GetIntSubfield(poCSAXField, "AXTY", 0);
     360          92 :             constexpr int AXIS_TYPE_GRAVITY_RELATED_DEPTH = 12;
     361          94 :             if (nAXTY != AXIS_TYPE_GRAVITY_RELATED_DEPTH &&
     362           2 :                 !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     363             :                     "Invalid value for AXTY field of CSAX idx %d.", nCSAXIdx)))
     364             :             {
     365           1 :                 return false;
     366             :             }
     367             : 
     368             :             //  Axis Unit of Measure
     369          91 :             const int nAXUM = poRecord->GetIntSubfield(poCSAXField, "AXUM", 0);
     370          91 :             constexpr int UOM_METRE = 4;
     371          93 :             if (nAXUM != UOM_METRE &&
     372           2 :                 !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     373             :                     "Invalid value for AXUM field of CSAX idx %d.", nCSAXIdx)))
     374             :             {
     375           1 :                 return false;
     376             :             }
     377             :         }
     378             : 
     379             :         // Read associated VDAT field
     380          95 :         CPLAssert(static_cast<size_t>(iCRSHIdx) < anVDATInstanceForCRSH.size());
     381          95 :         const int nVDATIdx = anVDATInstanceForCRSH[iCRSHIdx];
     382          95 :         if (nVDATIdx < 0 && !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     383             :                                 "No VDAT field for CRSH idx %d.", iCRSHIdx)))
     384             :         {
     385           1 :             return false;
     386             :         }
     387             : 
     388          94 :         const char *pszDTNM = nullptr;
     389          94 :         if (nVDATIdx >= 0)
     390             :         {
     391          85 :             const auto poVDATField = poRecord->FindField(VDAT_FIELD, nVDATIdx);
     392             : 
     393             :             // Vertical Datum Name
     394          85 :             pszDTNM = poRecord->GetStringSubfield(poVDATField, "DTNM", 0);
     395          85 :             if (!pszDTNM &&
     396           0 :                 !EMIT_ERROR(CPLSPrintf(
     397             :                     "No value for DTNM field of VDAT idx %d.", nVDATIdx)))
     398           0 :                 return false;
     399             : 
     400             :             // Vertical Datum Identifier
     401             :             const char *pszDTID =
     402          85 :                 poRecord->GetStringSubfield(poVDATField, "DTID", 0);
     403          85 :             if (!pszDTID &&
     404           0 :                 !EMIT_ERROR(CPLSPrintf(
     405             :                     "No value for DTID field of VDAT idx %d.", nVDATIdx)))
     406           0 :                 return false;
     407             : 
     408             :             // Datum Source
     409          85 :             const int nDTSR = poRecord->GetIntSubfield(poVDATField, "DTSR", 0);
     410          85 :             constexpr int SOURCE_FEATURE_CATALG = 2;
     411          87 :             if (nDTSR != SOURCE_FEATURE_CATALG &&
     412           2 :                 !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     413             :                     "Invalid value for DTSR field of VDAT idx %d.", nVDATIdx)))
     414           1 :                 return false;
     415             :         }
     416          93 :         if (!pszDTNM)
     417           9 :             pszDTNM = "(null)";
     418             : 
     419             :         // WKT CRS constraints
     420          93 :         if (strchr(pszCRNM, '"'))
     421             :         {
     422           2 :             if (!EMIT_ERROR_OR_WARNING(
     423             :                     "Double quote not allowed in vertical CRS name."))
     424           1 :                 return false;
     425           1 :             pszCRNM = "unknown";
     426             :         }
     427          92 :         if (strchr(pszDTNM, '"'))
     428             :         {
     429           2 :             if (!EMIT_ERROR_OR_WARNING(
     430             :                     "Double quote not allowed in vertical datum name."))
     431           1 :                 return false;
     432           1 :             pszDTNM = "unknown";
     433             :         }
     434             : 
     435          91 :         OGRSpatialReference oSRSVertical;
     436             :         const std::string osWKT = CPLSPrintf(
     437             :             "VERTCRS[\"%s depth\",VDATUM[\"%s\"],CS[vertical,1],"
     438             :             "AXIS[\"gravity-related depth (D)\",down,LENGTHUNIT[\"metre\",1]]]",
     439          91 :             pszCRNM, pszDTNM);
     440          91 :         if (oSRSVertical.importFromWkt(osWKT.c_str()) != OGRERR_NONE)
     441             :         {
     442           0 :             if (m_bStrict)
     443           0 :                 return false;
     444           0 :             continue;
     445             :         }
     446             : 
     447          91 :         OGRSpatialReference oCompoundCRS;
     448          91 :         oCompoundCRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     449         273 :         if (oCompoundCRS.SetCompoundCS(std::string(oSRSHorizontal.GetName())
     450          91 :                                            .append(" + ")
     451          91 :                                            .append(oSRSVertical.GetName())
     452             :                                            .c_str(),
     453             :                                        &oSRSHorizontal,
     454          91 :                                        &oSRSVertical) == OGRERR_NONE)
     455             :         {
     456          91 :             m_oMapSRS[nCRIX] = std::move(oCompoundCRS);
     457             :         }
     458           0 :         else if (m_bStrict)
     459             :         {
     460           0 :             return false;
     461             :         }
     462             :     }
     463             : 
     464         421 :     return true;
     465             : }