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