Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL Core
       4             :  * Purpose:  Read metadata from Pleiades imagery.
       5             :  * Author:   Alexander Lisovenko
       6             :  * Author:   Dmitry Baryshnikov, polimax@mail.ru
       7             :  *
       8             :  ******************************************************************************
       9             :  * Copyright (c) 2014-2018 NextGIS <info@nextgis.ru>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "reader_pleiades.h"
      16             : 
      17             : #include <cstddef>
      18             : #include <cstdio>
      19             : #include <cstring>
      20             : #include <ctime>
      21             : 
      22             : #include <string>
      23             : 
      24             : #include "cpl_conv.h"
      25             : #include "cpl_error.h"
      26             : #include "cpl_minixml.h"
      27             : #include "cpl_string.h"
      28             : #include "cpl_time.h"
      29             : 
      30             : /**
      31             :  * GDALMDReaderPleiades()
      32             :  */
      33       10184 : GDALMDReaderPleiades::GDALMDReaderPleiades(const char *pszPath,
      34       10184 :                                            char **papszSiblingFiles)
      35             :     : GDALMDReaderBase(pszPath, papszSiblingFiles), m_osBaseFilename(pszPath),
      36       10184 :       m_osIMDSourceFilename(CPLString()), m_osRPBSourceFilename(CPLString())
      37             : {
      38       10184 :     const CPLString osBaseName = CPLGetBasename(pszPath);
      39       10184 :     const size_t nBaseNameLen = osBaseName.size();
      40       10184 :     if (nBaseNameLen < 4 || nBaseNameLen > 511)
      41         768 :         return;
      42             : 
      43        9416 :     const CPLString osDirName = CPLGetDirname(pszPath);
      44             : 
      45             :     std::string osIMDSourceFilename = CPLFormFilename(
      46        9416 :         osDirName, CPLSPrintf("DIM_%s", osBaseName.c_str() + 4), "XML");
      47             :     std::string osRPBSourceFilename = CPLFormFilename(
      48        9416 :         osDirName, CPLSPrintf("RPC_%s", osBaseName.c_str() + 4), "XML");
      49             : 
      50             :     // find last underline
      51             :     char sBaseName[512];
      52        9416 :     size_t nLastUnderline = 0;
      53       76029 :     for (size_t i = 4; i < nBaseNameLen; i++)
      54             :     {
      55       66613 :         sBaseName[i - 4] = osBaseName[i];
      56       66613 :         if (osBaseName[i] == '_')
      57       10599 :             nLastUnderline = i - 4U;
      58             :     }
      59             : 
      60        9416 :     sBaseName[nLastUnderline] = 0;
      61             : 
      62             :     // Check if last 4 characters are fit in mask RjCj
      63             :     unsigned int iRow, iCol;
      64        9416 :     bool bHasRowColPart = nBaseNameLen > nLastUnderline + 5U;
      65        9416 :     if (!bHasRowColPart || sscanf(osBaseName.c_str() + nLastUnderline + 5U,
      66             :                                   "R%uC%u", &iRow, &iCol) != 2)
      67             :     {
      68        9402 :         return;
      69             :     }
      70             : 
      71             :     // Strip of suffix from PNEO products
      72          14 :     char *pszLastUnderScore = strrchr(sBaseName, '_');
      73          14 :     if (pszLastUnderScore &&
      74           3 :         (EQUAL(pszLastUnderScore, "_P") || EQUAL(pszLastUnderScore, "_RGB") ||
      75           3 :          EQUAL(pszLastUnderScore, "_NED")))
      76             :     {
      77           0 :         *pszLastUnderScore = 0;
      78             :     }
      79             : 
      80          14 :     if (CPLCheckForFile(&osIMDSourceFilename[0], papszSiblingFiles))
      81             :     {
      82           0 :         m_osIMDSourceFilename = std::move(osIMDSourceFilename);
      83             :     }
      84             :     else
      85             :     {
      86             :         osIMDSourceFilename =
      87          14 :             CPLFormFilename(osDirName, CPLSPrintf("DIM_%s", sBaseName), "XML");
      88          14 :         if (CPLCheckForFile(&osIMDSourceFilename[0], papszSiblingFiles))
      89             :         {
      90          14 :             m_osIMDSourceFilename = std::move(osIMDSourceFilename);
      91             :         }
      92             :     }
      93             : 
      94          14 :     if (CPLCheckForFile(&osRPBSourceFilename[0], papszSiblingFiles))
      95             :     {
      96           0 :         m_osRPBSourceFilename = std::move(osRPBSourceFilename);
      97             :     }
      98             :     else
      99             :     {
     100             :         osRPBSourceFilename =
     101          14 :             CPLFormFilename(osDirName, CPLSPrintf("RPC_%s", sBaseName), "XML");
     102          14 :         if (CPLCheckForFile(&osRPBSourceFilename[0], papszSiblingFiles))
     103             :         {
     104           4 :             m_osRPBSourceFilename = std::move(osRPBSourceFilename);
     105             :         }
     106             :     }
     107             : 
     108          14 :     if (!m_osIMDSourceFilename.empty())
     109          14 :         CPLDebug("MDReaderPleiades", "IMD Filename: %s",
     110             :                  m_osIMDSourceFilename.c_str());
     111          14 :     if (!m_osRPBSourceFilename.empty())
     112           4 :         CPLDebug("MDReaderPleiades", "RPB Filename: %s",
     113             :                  m_osRPBSourceFilename.c_str());
     114             : }
     115             : 
     116           6 : GDALMDReaderPleiades::GDALMDReaderPleiades()
     117           6 :     : GDALMDReaderBase(nullptr, nullptr)
     118             : {
     119           6 : }
     120             : 
     121             : GDALMDReaderPleiades *
     122           6 : GDALMDReaderPleiades::CreateReaderForRPC(const char *pszRPCSourceFilename)
     123             : {
     124           6 :     GDALMDReaderPleiades *poReader = new GDALMDReaderPleiades();
     125           6 :     poReader->m_osRPBSourceFilename = pszRPCSourceFilename;
     126           6 :     return poReader;
     127             : }
     128             : 
     129             : /**
     130             :  * ~GDALMDReaderPleiades()
     131             :  */
     132       15295 : GDALMDReaderPleiades::~GDALMDReaderPleiades()
     133             : {
     134       15295 : }
     135             : 
     136             : /**
     137             :  * HasRequiredFiles()
     138             :  */
     139       10184 : bool GDALMDReaderPleiades::HasRequiredFiles() const
     140             : {
     141       10184 :     if (!m_osIMDSourceFilename.empty())
     142          15 :         return true;
     143       10169 :     if (!m_osRPBSourceFilename.empty())
     144           0 :         return true;
     145             : 
     146       10169 :     return false;
     147             : }
     148             : 
     149             : /**
     150             :  * GetMetadataFiles()
     151             :  */
     152          15 : char **GDALMDReaderPleiades::GetMetadataFiles() const
     153             : {
     154          15 :     char **papszFileList = nullptr;
     155          15 :     if (!m_osIMDSourceFilename.empty())
     156          15 :         papszFileList = CSLAddString(papszFileList, m_osIMDSourceFilename);
     157          15 :     if (!m_osRPBSourceFilename.empty())
     158           4 :         papszFileList = CSLAddString(papszFileList, m_osRPBSourceFilename);
     159             : 
     160          15 :     return papszFileList;
     161             : }
     162             : 
     163             : /**
     164             :  * LoadMetadata()
     165             :  */
     166          25 : void GDALMDReaderPleiades::LoadMetadata()
     167             : {
     168          25 :     if (m_bIsMetadataLoad)
     169          11 :         return;
     170             : 
     171          14 :     if (!m_osIMDSourceFilename.empty())
     172             :     {
     173          14 :         CPLXMLNode *psNode = CPLParseXMLFile(m_osIMDSourceFilename);
     174             : 
     175          14 :         if (psNode != nullptr)
     176             :         {
     177          14 :             CPLXMLNode *psisdNode = CPLSearchXMLNode(psNode, "=Dimap_Document");
     178             : 
     179          14 :             if (psisdNode != nullptr)
     180             :             {
     181          14 :                 m_papszIMDMD = ReadXMLToList(psisdNode->psChild, m_papszIMDMD);
     182             :             }
     183          14 :             CPLDestroyXMLNode(psNode);
     184             :         }
     185             :     }
     186             : 
     187          14 :     if (!m_osRPBSourceFilename.empty())
     188             :     {
     189           4 :         m_papszRPCMD = LoadRPCXmlFile();
     190             :     }
     191             : 
     192          14 :     m_papszDEFAULTMD =
     193          14 :         CSLAddNameValue(m_papszDEFAULTMD, MD_NAME_MDTYPE, "DIMAP");
     194             : 
     195          14 :     m_bIsMetadataLoad = true;
     196             : 
     197          14 :     if (nullptr == m_papszIMDMD)
     198             :     {
     199           0 :         return;
     200             :     }
     201             : 
     202             :     // extract imagery metadata
     203          14 :     int nCounter = -1;
     204          28 :     const char *pszSatId1 = CSLFetchNameValue(
     205          14 :         m_papszIMDMD,
     206             :         "Dataset_Sources.Source_Identification.Strip_Source.MISSION");
     207          14 :     if (nullptr == pszSatId1)
     208             :     {
     209           0 :         nCounter = 1;
     210           0 :         for (int i = 0; i < 5; i++)
     211             :         {
     212           0 :             pszSatId1 = CSLFetchNameValue(
     213           0 :                 m_papszIMDMD,
     214             :                 CPLSPrintf("Dataset_Sources.Source_Identification_%d.Strip_"
     215             :                            "Source.MISSION",
     216             :                            nCounter));
     217           0 :             if (nullptr != pszSatId1)
     218           0 :                 break;
     219           0 :             nCounter++;
     220             :         }
     221             :     }
     222             : 
     223             :     const char *pszSatId2;
     224          14 :     if (nCounter == -1)
     225          14 :         pszSatId2 = CSLFetchNameValue(
     226          14 :             m_papszIMDMD,
     227             :             "Dataset_Sources.Source_Identification.Strip_Source.MISSION_INDEX");
     228             :     else
     229           0 :         pszSatId2 = CSLFetchNameValue(
     230           0 :             m_papszIMDMD, CPLSPrintf("Dataset_Sources.Source_Identification_%d."
     231             :                                      "Strip_Source.MISSION_INDEX",
     232             :                                      nCounter));
     233             : 
     234          14 :     if (nullptr != pszSatId1 && nullptr != pszSatId2)
     235             :     {
     236          28 :         m_papszIMAGERYMD = CSLAddNameValue(
     237             :             m_papszIMAGERYMD, MD_NAME_SATELLITE,
     238          28 :             CPLSPrintf("%s %s", CPLStripQuotes(pszSatId1).c_str(),
     239          28 :                        CPLStripQuotes(pszSatId2).c_str()));
     240             :     }
     241           0 :     else if (nullptr != pszSatId1 && nullptr == pszSatId2)
     242             :     {
     243           0 :         m_papszIMAGERYMD = CSLAddNameValue(m_papszIMAGERYMD, MD_NAME_SATELLITE,
     244           0 :                                            CPLStripQuotes(pszSatId1));
     245             :     }
     246           0 :     else if (nullptr == pszSatId1 && nullptr != pszSatId2)
     247             :     {
     248           0 :         m_papszIMAGERYMD = CSLAddNameValue(m_papszIMAGERYMD, MD_NAME_SATELLITE,
     249           0 :                                            CPLStripQuotes(pszSatId2));
     250             :     }
     251             : 
     252             :     const char *pszDate;
     253          14 :     if (nCounter == -1)
     254          14 :         pszDate = CSLFetchNameValue(
     255          14 :             m_papszIMDMD,
     256             :             "Dataset_Sources.Source_Identification.Strip_Source.IMAGING_DATE");
     257             :     else
     258           0 :         pszDate = CSLFetchNameValue(
     259           0 :             m_papszIMDMD, CPLSPrintf("Dataset_Sources.Source_Identification_%d."
     260             :                                      "Strip_Source.IMAGING_DATE",
     261             :                                      nCounter));
     262             : 
     263          14 :     if (nullptr != pszDate)
     264             :     {
     265             :         const char *pszTime;
     266          14 :         if (nCounter == -1)
     267          14 :             pszTime = CSLFetchNameValue(m_papszIMDMD,
     268             :                                         "Dataset_Sources.Source_Identification."
     269             :                                         "Strip_Source.IMAGING_TIME");
     270             :         else
     271           0 :             pszTime = CSLFetchNameValue(
     272           0 :                 m_papszIMDMD,
     273             :                 CPLSPrintf("Dataset_Sources.Source_Identification_%d.Strip_"
     274             :                            "Source.IMAGING_TIME",
     275             :                            nCounter));
     276             : 
     277          14 :         if (nullptr == pszTime)
     278           0 :             pszTime = "00:00:00.0Z";
     279             : 
     280             :         char buffer[80];
     281             :         GIntBig timeMid =
     282          14 :             GetAcquisitionTimeFromString(CPLSPrintf("%sT%s", pszDate, pszTime));
     283             :         struct tm tmBuf;
     284          14 :         strftime(buffer, 80, MD_DATETIMEFORMAT,
     285          14 :                  CPLUnixTimeToYMDHMS(timeMid, &tmBuf));
     286          14 :         m_papszIMAGERYMD =
     287          14 :             CSLAddNameValue(m_papszIMAGERYMD, MD_NAME_ACQDATETIME, buffer);
     288             :     }
     289             : 
     290          14 :     m_papszIMAGERYMD =
     291          14 :         CSLAddNameValue(m_papszIMAGERYMD, MD_NAME_CLOUDCOVER, MD_CLOUDCOVER_NA);
     292             : }
     293             : 
     294             : /**
     295             :  * LoadRPCXmlFile()
     296             :  */
     297             : 
     298             : static const char *const apszRPBMapPleiades[] = {
     299             :     RPC_LINE_OFF,     "RFM_Validity.LINE_OFF",  // do not change order !
     300             :     RPC_SAMP_OFF,     "RFM_Validity.SAMP_OFF",  // do not change order !
     301             :     RPC_LAT_OFF,      "RFM_Validity.LAT_OFF",
     302             :     RPC_LONG_OFF,     "RFM_Validity.LONG_OFF",
     303             :     RPC_HEIGHT_OFF,   "RFM_Validity.HEIGHT_OFF",
     304             :     RPC_LINE_SCALE,   "RFM_Validity.LINE_SCALE",
     305             :     RPC_SAMP_SCALE,   "RFM_Validity.SAMP_SCALE",
     306             :     RPC_LAT_SCALE,    "RFM_Validity.LAT_SCALE",
     307             :     RPC_LONG_SCALE,   "RFM_Validity.LONG_SCALE",
     308             :     RPC_HEIGHT_SCALE, "RFM_Validity.HEIGHT_SCALE",
     309             :     nullptr,          nullptr};
     310             : 
     311             : static const char *const apszRPCTXT20ValItemsPleiades[] = {
     312             :     RPC_LINE_NUM_COEFF, RPC_LINE_DEN_COEFF, RPC_SAMP_NUM_COEFF,
     313             :     RPC_SAMP_DEN_COEFF, nullptr};
     314             : 
     315          10 : char **GDALMDReaderPleiades::LoadRPCXmlFile()
     316             : {
     317          10 :     CPLXMLNode *pNode = CPLParseXMLFile(m_osRPBSourceFilename);
     318             : 
     319          10 :     if (nullptr == pNode)
     320           0 :         return nullptr;
     321             : 
     322             :     // search Global_RFM
     323          10 :     char **papszRawRPCList = nullptr;
     324          10 :     CPLXMLNode *pGRFMNode = CPLSearchXMLNode(pNode, "=Global_RFM");
     325             : 
     326          10 :     if (pGRFMNode != nullptr)
     327             :     {
     328          10 :         papszRawRPCList = ReadXMLToList(pGRFMNode->psChild, papszRawRPCList);
     329             :     }
     330             :     else
     331             :     {
     332           0 :         pGRFMNode = CPLSearchXMLNode(pNode, "=Rational_Function_Model");
     333             : 
     334           0 :         if (pGRFMNode != nullptr)
     335             :         {
     336             :             papszRawRPCList =
     337           0 :                 ReadXMLToList(pGRFMNode->psChild, papszRawRPCList);
     338             :         }
     339             :     }
     340             : 
     341          10 :     if (nullptr == papszRawRPCList)
     342             :     {
     343           0 :         CPLDestroyXMLNode(pNode);
     344           0 :         return nullptr;
     345             :     }
     346             : 
     347             :     // If we are not the top-left tile, then we must shift LINE_OFF and SAMP_OFF
     348          10 :     int nLineOffShift = 0;
     349          10 :     int nPixelOffShift = 0;
     350          10 :     for (int i = 1; TRUE; i++)
     351             :     {
     352          11 :         CPLString osKey;
     353             :         osKey.Printf("Raster_Data.Data_Access.Data_Files.Data_File_%d.DATA_"
     354             :                      "FILE_PATH.href",
     355          11 :                      i);
     356          11 :         const char *pszHref = CSLFetchNameValue(m_papszIMDMD, osKey);
     357          11 :         if (pszHref == nullptr)
     358           9 :             break;
     359           2 :         if (strcmp(CPLGetFilename(pszHref), CPLGetFilename(m_osBaseFilename)) ==
     360             :             0)
     361             :         {
     362             :             osKey.Printf(
     363           1 :                 "Raster_Data.Data_Access.Data_Files.Data_File_%d.tile_C", i);
     364           1 :             const char *pszC = CSLFetchNameValue(m_papszIMDMD, osKey);
     365             :             osKey.Printf(
     366           1 :                 "Raster_Data.Data_Access.Data_Files.Data_File_%d.tile_R", i);
     367           1 :             const char *pszR = CSLFetchNameValue(m_papszIMDMD, osKey);
     368           2 :             const char *pszTileWidth = CSLFetchNameValue(
     369           1 :                 m_papszIMDMD, "Raster_Data.Raster_Dimensions.Tile_Set.Regular_"
     370             :                               "Tiling.NTILES_SIZE.ncols");
     371           2 :             const char *pszTileHeight = CSLFetchNameValue(
     372           1 :                 m_papszIMDMD, "Raster_Data.Raster_Dimensions.Tile_Set.Regular_"
     373             :                               "Tiling.NTILES_SIZE.nrows");
     374             :             const char *pszOVERLAP_COL =
     375           1 :                 CSLFetchNameValueDef(m_papszIMDMD,
     376             :                                      "Raster_Data.Raster_Dimensions.Tile_Set."
     377             :                                      "Regular_Tiling.OVERLAP_COL",
     378             :                                      "0");
     379             :             const char *pszOVERLAP_ROW =
     380           1 :                 CSLFetchNameValueDef(m_papszIMDMD,
     381             :                                      "Raster_Data.Raster_Dimensions.Tile_Set."
     382             :                                      "Regular_Tiling.OVERLAP_ROW",
     383             :                                      "0");
     384             : 
     385           1 :             if (pszC && pszR && pszTileWidth && pszTileHeight &&
     386           1 :                 atoi(pszOVERLAP_COL) == 0 && atoi(pszOVERLAP_ROW) == 0)
     387             :             {
     388           1 :                 nLineOffShift = -(atoi(pszR) - 1) * atoi(pszTileHeight);
     389           1 :                 nPixelOffShift = -(atoi(pszC) - 1) * atoi(pszTileWidth);
     390             :             }
     391           1 :             break;
     392             :         }
     393           1 :     }
     394             : 
     395             :     // SPOT and PHR sensors use 1,1 as their upper left corner pixel convention
     396             :     // for RPCs which is non standard. This was fixed with PNEO which correctly
     397             :     // assumes 0,0.
     398             :     // Precompute the offset that will be applied to LINE_OFF and SAMP_OFF
     399             :     // in order to use the RPCs with the standard 0,0 convention
     400             :     double topleftOffset;
     401          10 :     CPLXMLNode *psDoc = CPLGetXMLNode(pNode, "=Dimap_Document");
     402          10 :     if (!psDoc)
     403           0 :         psDoc = CPLGetXMLNode(pNode, "=PHR_DIMAP_Document");
     404          10 :     const char *pszMetadataProfile = CPLGetXMLValue(
     405             :         psDoc, "Metadata_Identification.METADATA_PROFILE", "PHR_SENSOR");
     406          10 :     if (EQUAL(pszMetadataProfile, "PHR_SENSOR") ||
     407           2 :         EQUAL(pszMetadataProfile, "S7_SENSOR") ||
     408           2 :         EQUAL(pszMetadataProfile, "S6_SENSOR"))
     409             :     {
     410           8 :         topleftOffset = 1;
     411             :     }
     412           2 :     else if (EQUAL(pszMetadataProfile, "PNEO_SENSOR"))
     413             :     {
     414           2 :         topleftOffset = 0;
     415             :     }
     416             :     else
     417             :     {
     418             :         //CPLError(CE_Warning, CPLE_AppDefined,
     419             :         //         "Unknown RPC Metadata Profile: %s. Assuming PHR_SENSOR",
     420             :         //         pszMetadataProfile);
     421           0 :         topleftOffset = 1;
     422             :     }
     423             : 
     424             :     // format list
     425          10 :     char **papszRPB = nullptr;
     426         110 :     for (int i = 0; apszRPBMapPleiades[i] != nullptr; i += 2)
     427             :     {
     428             :         const char *pszValue =
     429         100 :             CSLFetchNameValue(papszRawRPCList, apszRPBMapPleiades[i + 1]);
     430         100 :         if ((i == 0 || i == 2) && pszValue)  //i.e. LINE_OFF or SAMP_OFF
     431             :         {
     432          20 :             CPLString osField;
     433          20 :             double dfVal = CPLAtofM(pszValue) - topleftOffset;
     434          20 :             if (i == 0)
     435          10 :                 dfVal += nLineOffShift;
     436             :             else
     437          10 :                 dfVal += nPixelOffShift;
     438          20 :             osField.Printf("%.15g", dfVal);
     439             :             papszRPB =
     440          20 :                 CSLAddNameValue(papszRPB, apszRPBMapPleiades[i], osField);
     441             :         }
     442             :         else
     443             :         {
     444             :             papszRPB =
     445          80 :                 CSLAddNameValue(papszRPB, apszRPBMapPleiades[i], pszValue);
     446             :         }
     447             :     }
     448             : 
     449             :     // merge coefficients
     450          50 :     for (int i = 0; apszRPCTXT20ValItemsPleiades[i] != nullptr; i++)
     451             :     {
     452          40 :         CPLString value;
     453         840 :         for (int j = 1; j < 21; j++)
     454             :         {
     455             :             // We want to use the Inverse_Model
     456             :             // Quoting PleiadesUserGuideV2-1012.pdf:
     457             :             // """When using the inverse model (ground --> image), the user
     458             :             // supplies geographic coordinates (lon, lat) and an altitude
     459             :             // (alt)"""
     460         800 :             const char *pszValue = CSLFetchNameValue(
     461             :                 papszRawRPCList,
     462             :                 CPLSPrintf("Inverse_Model.%s_%d",
     463         800 :                            apszRPCTXT20ValItemsPleiades[i], j));
     464         800 :             if (nullptr != pszValue)
     465             :             {
     466         640 :                 value = value + " " + CPLString(pszValue);
     467             :             }
     468             :             else
     469             :             {
     470         160 :                 pszValue = CSLFetchNameValue(
     471             :                     papszRawRPCList,
     472             :                     CPLSPrintf("GroundtoImage_Values.%s_%d",
     473         160 :                                apszRPCTXT20ValItemsPleiades[i], j));
     474         160 :                 if (nullptr != pszValue)
     475             :                 {
     476         160 :                     value = value + " " + CPLString(pszValue);
     477             :                 }
     478             :             }
     479             :         }
     480             :         papszRPB =
     481          40 :             CSLAddNameValue(papszRPB, apszRPCTXT20ValItemsPleiades[i], value);
     482             :     }
     483             : 
     484          10 :     CSLDestroy(papszRawRPCList);
     485          10 :     CPLDestroyXMLNode(pNode);
     486          10 :     return papszRPB;
     487             : }