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