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             : /*                         GetPointLayerName()                          */
      20             : /************************************************************************/
      21             : 
      22             : /* static */ std::string
      23         560 : OGRS101Reader::GetPointLayerName(const OGRSpatialReference &oSRS)
      24             : {
      25         560 :     return oSRS.GetAxesCount() == 2
      26             :                ? "Point2D"
      27        1680 :                : CPLSPrintf("Point3D_%s", LaunderCRSName(oSRS).c_str());
      28             : }
      29             : 
      30             : /************************************************************************/
      31             : /*                      CreatePointFeatureDefns()                       */
      32             : /************************************************************************/
      33             : 
      34             : /** Create the feature definition(s) for the Point layer(s)
      35             :  *
      36             :  * There is a layer per CRS used by points.
      37             :  */
      38         353 : bool OGRS101Reader::CreatePointFeatureDefns()
      39             : {
      40         353 :     bool bError = false;
      41         353 :     m_oMapCRSIdToPointRecordIdx = CreateMapCRSIdToRecordIdxForPoints(bError);
      42         353 :     if (bError)
      43           2 :         return false;
      44         575 :     for (const auto &[nCRSId, anRecordIdx] : m_oMapCRSIdToPointRecordIdx)
      45             :     {
      46         228 :         const auto &oSRS = m_oMapSRS[nCRSId];
      47         228 :         const bool bIs2D = nCRSId == HORIZONTAL_CRS_ID;
      48             :         auto poFDefn = OGRFeatureDefnRefCountedPtr::makeInstance(
      49         228 :             GetPointLayerName(oSRS).c_str());
      50         228 :         poFDefn->SetGeomType(bIs2D ? wkbPoint : wkbPoint25D);
      51         456 :         poFDefn->GetGeomFieldDefn(0)->SetSpatialRef(
      52         456 :             OGRSpatialReferenceRefCountedPtr::makeClone(&oSRS).get());
      53         228 :         poFDefn->GetGeomFieldDefn(0)->SetCoordinatePrecision(
      54         228 :             m_coordinatePrecision);
      55             :         {
      56         456 :             OGRFieldDefn oFieldDefn(OGR_FIELD_NAME_RECORD_ID, OFTInteger);
      57         228 :             poFDefn->AddFieldDefn(&oFieldDefn);
      58             :         }
      59             :         {
      60         456 :             OGRFieldDefn oFieldDefn(OGR_FIELD_NAME_RECORD_VERSION, OFTInteger);
      61         228 :             poFDefn->AddFieldDefn(&oFieldDefn);
      62             :         }
      63             : 
      64         228 :         if (!InferFeatureDefn(m_oPointRecordIndex, PRID_FIELD, INAS_FIELD,
      65         228 :                               anRecordIdx, *poFDefn, m_oMapFieldDomains))
      66             :         {
      67           4 :             return false;
      68             :         }
      69         224 :         m_oMapPointFeatureDefn[nCRSId] = std::move(poFDefn);
      70             :     }
      71             : 
      72         347 :     return true;
      73             : }
      74             : 
      75             : /************************************************************************/
      76             : /*                       GetCRSIdForPointRecord()                       */
      77             : /************************************************************************/
      78             : 
      79             : /** Return the CRS id for a given Point record, or INVALID_CRS_ID on error */
      80             : OGRS101Reader::CRSId
      81         523 : OGRS101Reader::GetCRSIdForPointRecord(const DDFRecord *poRecord, int iRecord,
      82             :                                       int nRecordID) const
      83             : {
      84         523 :     if (nRecordID < 0)
      85         427 :         nRecordID = poRecord->GetIntSubfield(PRID_FIELD, 0, RCID_SUBFIELD, 0);
      86             : 
      87          14 :     const auto GetErrorContext = [iRecord, nRecordID]()
      88             :     {
      89           7 :         if (iRecord >= 0)
      90           7 :             return CPLSPrintf("Record index=%d of PRID", iRecord);
      91             :         else
      92           0 :             return CPLSPrintf("Record ID=%d of PRID", nRecordID);
      93         523 :     };
      94             : 
      95         523 :     if (poRecord->FindField(C3IT_FIELD))
      96             :     {
      97             :         const CRSId nVCID =
      98          43 :             poRecord->GetIntSubfield(C3IT_FIELD, 0, VCID_SUBFIELD, 0);
      99          43 :         if (nVCID == HORIZONTAL_CRS_ID)
     100             :         {
     101           2 :             CPL_IGNORE_RET_VAL(EMIT_ERROR_OR_WARNING(
     102             :                 CPLSPrintf("%s: VCID subfield = %d of C3IT "
     103             :                            "field points to a non-3D CRS.",
     104             :                            GetErrorContext(), static_cast<int>(nVCID))));
     105           2 :             return INVALID_CRS_ID;
     106             :         }
     107          41 :         else if (!cpl::contains(m_oMapSRS, nVCID))
     108             :         {
     109           2 :             CPL_IGNORE_RET_VAL(EMIT_ERROR_OR_WARNING(
     110             :                 CPLSPrintf("%s: Unknown value %d for VCID subfield of C3IT "
     111             :                            "field.",
     112             :                            GetErrorContext(), static_cast<int>(nVCID))));
     113           2 :             return INVALID_CRS_ID;
     114             :         }
     115             :         else
     116             :         {
     117          39 :             return nVCID;
     118             :         }
     119             :     }
     120         480 :     else if (poRecord->FindField(C2IT_FIELD))
     121             :     {
     122         477 :         return HORIZONTAL_CRS_ID;
     123             :     }
     124             :     else
     125             :     {
     126           3 :         CPL_IGNORE_RET_VAL(EMIT_ERROR_OR_WARNING(
     127             :             CPLSPrintf("%s: No C2IT or C3IT field found.", GetErrorContext())));
     128           3 :         return INVALID_CRS_ID;
     129             :     }
     130             : }
     131             : 
     132             : /************************************************************************/
     133             : /*                 CreateMapCRSIdToRecordIdxForPoints()                 */
     134             : /************************************************************************/
     135             : 
     136             : /** Browse through m_oPointRecordIndex to identify which record belongs to
     137             :  * each CRS and create a map from each CRS id to the record indices that use
     138             :  * it.
     139             :  */
     140             : std::map<OGRS101Reader::CRSId, std::vector<int>>
     141         353 : OGRS101Reader::CreateMapCRSIdToRecordIdxForPoints(bool &bError) const
     142             : {
     143         706 :     std::map<OGRS101Reader::CRSId, std::vector<int>> map;
     144             : 
     145         353 :     const int nRecords = m_oPointRecordIndex.GetCount();
     146         778 :     for (int iRecord = 0; iRecord < nRecords; ++iRecord)
     147             :     {
     148         427 :         const auto poRecord = m_oPointRecordIndex.GetByIndex(iRecord);
     149             :         const CRSId nCRSId =
     150         427 :             GetCRSIdForPointRecord(poRecord, iRecord, /* nRecordID = */ -1);
     151         427 :         if (nCRSId == INVALID_CRS_ID)
     152             :         {
     153           7 :             if (m_bStrict)
     154             :             {
     155           2 :                 bError = true;
     156           2 :                 return {};
     157             :             }
     158             :         }
     159             :         else
     160             :         {
     161         420 :             map[nCRSId].push_back(iRecord);
     162             :         }
     163             :     }
     164             : 
     165         351 :     return map;
     166             : }
     167             : 
     168             : /************************************************************************/
     169             : /*                         ReadPointGeometry()                          */
     170             : /************************************************************************/
     171             : 
     172       15724 : std::unique_ptr<OGRPoint> OGRS101Reader::ReadPointGeometryInternal(
     173             :     const DDFRecord *poRecord, int iRecord, int nRecordID, int iPnt,
     174             :     const OGRSpatialReference *poSRS, const bool bIs3D,
     175             :     const DDFField *poCoordField, const char *pszRecordFieldName) const
     176             : {
     177          58 :     const auto GetErrorContext = [iRecord, nRecordID, pszRecordFieldName]()
     178             :     {
     179          29 :         if (iRecord >= 0)
     180          22 :             return CPLSPrintf("Record index=%d of %s", iRecord,
     181          22 :                               pszRecordFieldName);
     182             :         else
     183           7 :             return CPLSPrintf("Record ID=%d of %s", nRecordID,
     184           7 :                               pszRecordFieldName);
     185       15724 :     };
     186             : 
     187       15724 :     if (!poCoordField)
     188             :     {
     189           0 :         CPL_IGNORE_RET_VAL(EMIT_ERROR_OR_WARNING(CPLSPrintf(
     190             :             "%s: cannot find coordinate field.", GetErrorContext())));
     191           0 :         return nullptr;
     192             :     }
     193             : 
     194       15724 :     int bSuccess = false;
     195             :     const double dfX =
     196       15724 :         poRecord->GetIntSubfield(poCoordField, XCOO_SUBFIELD, iPnt, &bSuccess) /
     197       15724 :             static_cast<double>(m_nXScale) +
     198       15724 :         m_dfXShift;
     199       15724 :     if (!bSuccess &&
     200           0 :         !EMIT_ERROR_OR_WARNING(CPLSPrintf("%s: cannot read %s.",
     201             :                                           GetErrorContext(), XCOO_SUBFIELD)))
     202             :     {
     203           0 :         return nullptr;
     204             :     }
     205             :     const double dfY =
     206       15724 :         poRecord->GetIntSubfield(poCoordField, YCOO_SUBFIELD, iPnt, &bSuccess) /
     207       15724 :             static_cast<double>(m_nYScale) +
     208       15724 :         m_dfYShift;
     209       15724 :     if (!bSuccess &&
     210           0 :         !EMIT_ERROR_OR_WARNING(CPLSPrintf("%s: cannot read %s.",
     211             :                                           GetErrorContext(), YCOO_SUBFIELD)))
     212             :     {
     213           0 :         return nullptr;
     214             :     }
     215       15753 :     if (poSRS && !(std::fabs(dfX) <= 180 && std::fabs(dfY) <= 90) &&
     216          29 :         poSRS->IsGeographic())
     217             :     {
     218          29 :         if (!EMIT_ERROR_OR_WARNING(
     219             :                 CPLSPrintf("%s: wrong coordinate value: lon=%f, lat=%f.",
     220             :                            GetErrorContext(), dfX, dfY)))
     221             :         {
     222           8 :             return nullptr;
     223             :         }
     224             :     }
     225       15716 :     std::unique_ptr<OGRPoint> poPoint;
     226       15716 :     if (bIs3D)
     227             :     {
     228         390 :         const double dfZ = poRecord->GetIntSubfield(poCoordField, ZCOO_SUBFIELD,
     229         390 :                                                     iPnt, &bSuccess) /
     230         390 :                                static_cast<double>(m_nZScale) +
     231         390 :                            m_dfZShift;
     232         390 :         if (!bSuccess &&
     233           0 :             !EMIT_ERROR_OR_WARNING(CPLSPrintf(
     234             :                 "%s: cannot read %s.", GetErrorContext(), ZCOO_SUBFIELD)))
     235             :         {
     236           0 :             return nullptr;
     237             :         }
     238         390 :         poPoint = std::make_unique<OGRPoint>(dfX, dfY, dfZ);
     239             :     }
     240             :     else
     241             :     {
     242       15326 :         poPoint = std::make_unique<OGRPoint>(dfX, dfY);
     243             :     }
     244       15716 :     poPoint->assignSpatialReference(poSRS);
     245       15716 :     return poPoint;
     246             : }
     247             : 
     248             : std::unique_ptr<OGRPoint>
     249        1447 : OGRS101Reader::ReadPointGeometry(const DDFRecord *poRecord, int iRecord,
     250             :                                  int nRecordID,
     251             :                                  const OGRSpatialReference *poSRS) const
     252             : {
     253        1447 :     const DDFField *poCoordField = poRecord->FindField(C3IT_FIELD);
     254        1447 :     const bool bIs3D = poCoordField != nullptr;
     255        1447 :     if (!poCoordField)
     256        1268 :         poCoordField = poRecord->FindField(C2IT_FIELD);
     257             : 
     258             :     return ReadPointGeometryInternal(poRecord, iRecord, nRecordID,
     259             :                                      /* iPnt = */ 0, poSRS, bIs3D, poCoordField,
     260        1447 :                                      PRID_FIELD);
     261             : }
     262             : 
     263             : /************************************************************************/
     264             : /*                          FillFeaturePoint()                          */
     265             : /************************************************************************/
     266             : 
     267             : /** Fill the content of the provided feature from the identified record
     268             :  * (of m_oPointRecordIndex).
     269             :  */
     270        1146 : bool OGRS101Reader::FillFeaturePoint(const DDFRecordIndex &oIndex, int iRecord,
     271             :                                      OGRFeature &oFeature) const
     272             : {
     273        1146 :     const auto poRecord = oIndex.GetByIndex(iRecord);
     274        1146 :     CPLAssert(poRecord);
     275             : 
     276             :     const OGRSpatialReference *poSRS =
     277        1146 :         oFeature.GetDefnRef()->GetGeomFieldDefn(0)->GetSpatialRef();
     278             :     auto poPoint =
     279        2292 :         ReadPointGeometry(poRecord, iRecord, /* nRecordID = */ -1, poSRS);
     280        1146 :     if (poPoint)
     281        1141 :         oFeature.SetGeometry(std::move(poPoint));
     282           5 :     else if (m_bStrict)
     283           5 :         return false;
     284             : 
     285        2282 :     return FillFeatureAttributes(oIndex, iRecord, INAS_FIELD, oFeature) &&
     286        1141 :            FillFeatureWithNonAttrAssocSubfields(poRecord, iRecord, INAS_FIELD,
     287        1141 :                                                 oFeature);
     288             : }
     289             : 
     290             : /************************************************************************/
     291             : /*                      ProcessUpdateRecordPoint()                      */
     292             : /************************************************************************/
     293             : 
     294             : /** Updates the geometry part of poTargetRecord with poUpdateRecord */
     295          11 : bool OGRS101Reader::ProcessUpdateRecordPoint(const DDFRecord *poUpdateRecord,
     296             :                                              DDFRecord *poTargetRecord) const
     297             : {
     298          11 :     const auto poIDField = poUpdateRecord->GetField(0);
     299          11 :     CPLAssert(poIDField);
     300             : 
     301             :     // Record name
     302             :     const RecordName nRCNM =
     303          11 :         poUpdateRecord->GetIntSubfield(poIDField, RCNM_SUBFIELD, 0);
     304             : 
     305             :     // Record identifier
     306             :     const int nRCID =
     307          11 :         poUpdateRecord->GetIntSubfield(poIDField, RCID_SUBFIELD, 0);
     308             : 
     309          11 :     if (const auto poC2ITFieldUpdate = poUpdateRecord->FindField(C2IT_FIELD))
     310             :     {
     311           8 :         auto poC2ITFieldTarget = poTargetRecord->FindField(C2IT_FIELD);
     312           8 :         if (!poC2ITFieldTarget)
     313             :         {
     314           0 :             return EMIT_ERROR_OR_WARNING(CPLSPrintf(
     315             :                 "%s, RCNM=%d, RCID=%d: cannot find C2IT field in target "
     316             :                 "record",
     317             :                 m_osFilename.c_str(), static_cast<int>(nRCNM), nRCID));
     318             :         }
     319           8 :         if (*(poC2ITFieldTarget->GetFieldDefn()) !=
     320           8 :             *(poC2ITFieldUpdate->GetFieldDefn()))
     321             :         {
     322           0 :             return EMIT_ERROR("C2IT field definitions of update and target "
     323             :                               "records are different");
     324             :         }
     325          24 :         for (const char *pszSubFieldName : {XCOO_SUBFIELD, YCOO_SUBFIELD})
     326             :         {
     327          16 :             poTargetRecord->SetIntSubfield(
     328             :                 C2IT_FIELD, 0, pszSubFieldName, 0,
     329             :                 poUpdateRecord->GetIntSubfield(C2IT_FIELD, 0, pszSubFieldName,
     330             :                                                0));
     331             :         }
     332             :     }
     333           3 :     else if (const auto poC3ITFieldUpdate =
     334           3 :                  poUpdateRecord->FindField(C3IT_FIELD))
     335             :     {
     336           3 :         auto poC3ITFieldTarget = poTargetRecord->FindField(C3IT_FIELD);
     337           3 :         if (!poC3ITFieldTarget)
     338             :         {
     339           1 :             return EMIT_ERROR_OR_WARNING(CPLSPrintf(
     340             :                 "%s, RCNM=%d, RCID=%d: cannot find C3IT field in target "
     341             :                 "record",
     342             :                 m_osFilename.c_str(), static_cast<int>(nRCNM), nRCID));
     343             :         }
     344           2 :         if (*(poC3ITFieldTarget->GetFieldDefn()) !=
     345           2 :             *(poC3ITFieldUpdate->GetFieldDefn()))
     346             :         {
     347           0 :             return EMIT_ERROR("C3IT field definitions of update and target "
     348             :                               "records are different");
     349             :         }
     350           6 :         for (const char *pszSubFieldName :
     351           8 :              {XCOO_SUBFIELD, YCOO_SUBFIELD, ZCOO_SUBFIELD})
     352             :         {
     353           6 :             poTargetRecord->SetIntSubfield(
     354             :                 C3IT_FIELD, 0, pszSubFieldName, 0,
     355             :                 poUpdateRecord->GetIntSubfield(C3IT_FIELD, 0, pszSubFieldName,
     356             :                                                0));
     357             :         }
     358             :     }
     359             : 
     360          10 :     return true;
     361             : }