Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Hierarchical Data Format Release 5 (HDF5)
4 : * Purpose: Read S100 bathymetric datasets.
5 : * Author: Even Rouault <even dot rouault at spatialys dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #ifdef _POSIX_C_SOURCE
14 : #undef _POSIX_C_SOURCE
15 : #endif
16 :
17 : #include "cpl_time.h"
18 :
19 : #include "s100.h"
20 : #include "hdf5dataset.h"
21 : #include "gh5_convenience.h"
22 :
23 : #include "proj.h"
24 : #include "proj_experimental.h"
25 : #include "proj_constants.h"
26 : #include "ogr_proj_p.h"
27 :
28 : #include <algorithm>
29 : #include <array>
30 : #include <cmath>
31 : #include <limits>
32 :
33 : /************************************************************************/
34 : /* S100BaseDataset() */
35 : /************************************************************************/
36 :
37 224 : S100BaseDataset::S100BaseDataset(const std::string &osFilename)
38 224 : : m_osFilename(osFilename)
39 : {
40 224 : }
41 :
42 : /************************************************************************/
43 : /* Init() */
44 : /************************************************************************/
45 :
46 224 : bool S100BaseDataset::Init()
47 : {
48 : // Open the file as an HDF5 file.
49 224 : hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
50 224 : H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr);
51 224 : hid_t hHDF5 = H5Fopen(m_osFilename.c_str(), H5F_ACC_RDONLY, fapl);
52 224 : H5Pclose(fapl);
53 224 : if (hHDF5 < 0)
54 0 : return false;
55 :
56 448 : auto poSharedResources = GDAL::HDF5SharedResources::Create(m_osFilename);
57 224 : poSharedResources->m_hHDF5 = hHDF5;
58 :
59 224 : m_poRootGroup = HDF5Dataset::OpenGroup(poSharedResources);
60 224 : if (m_poRootGroup == nullptr)
61 0 : return false;
62 :
63 224 : S100ReadSRS(m_poRootGroup.get(), m_oSRS);
64 :
65 224 : S100ReadVerticalDatum(this, m_poRootGroup.get());
66 :
67 : m_osMetadataFile =
68 224 : S100ReadMetadata(this, m_osFilename, m_poRootGroup.get());
69 :
70 224 : return true;
71 : }
72 :
73 : /************************************************************************/
74 : /* GetGeoTransform() */
75 : /************************************************************************/
76 :
77 65 : CPLErr S100BaseDataset::GetGeoTransform(GDALGeoTransform >) const
78 :
79 : {
80 65 : if (m_bHasGT)
81 : {
82 65 : gt = m_gt;
83 65 : return CE_None;
84 : }
85 :
86 0 : return GDALPamDataset::GetGeoTransform(gt);
87 : }
88 :
89 : /************************************************************************/
90 : /* GetSpatialRef() */
91 : /************************************************************************/
92 :
93 97 : const OGRSpatialReference *S100BaseDataset::GetSpatialRef() const
94 : {
95 97 : if (!m_oSRS.IsEmpty())
96 97 : return &m_oSRS;
97 0 : return GDALPamDataset::GetSpatialRef();
98 : }
99 :
100 : /************************************************************************/
101 : /* GetFileList() */
102 : /************************************************************************/
103 :
104 12 : char **S100BaseDataset::GetFileList()
105 : {
106 12 : char **papszFileList = GDALPamDataset::GetFileList();
107 12 : if (!m_osMetadataFile.empty())
108 4 : papszFileList = CSLAddString(papszFileList, m_osMetadataFile.c_str());
109 12 : return papszFileList;
110 : }
111 :
112 : /************************************************************************/
113 : /* S100ReadSRS() */
114 : /************************************************************************/
115 :
116 : constexpr int PROJECTION_METHOD_MERCATOR = 9805;
117 : static_assert(PROJECTION_METHOD_MERCATOR ==
118 : EPSG_CODE_METHOD_MERCATOR_VARIANT_B);
119 : constexpr int PROJECTION_METHOD_TRANSVERSE_MERCATOR = 9807;
120 : static_assert(PROJECTION_METHOD_TRANSVERSE_MERCATOR ==
121 : EPSG_CODE_METHOD_TRANSVERSE_MERCATOR);
122 : constexpr int PROJECTION_METHOD_OBLIQUE_MERCATOR = 9815;
123 : static_assert(PROJECTION_METHOD_OBLIQUE_MERCATOR ==
124 : EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_B);
125 : constexpr int PROJECTION_METHOD_HOTINE_OBLIQUE_MERCATOR = 9812;
126 : static_assert(PROJECTION_METHOD_HOTINE_OBLIQUE_MERCATOR ==
127 : EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_A);
128 : constexpr int PROJECTION_METHOD_LCC_1SP = 9801;
129 : static_assert(PROJECTION_METHOD_LCC_1SP ==
130 : EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_1SP);
131 : constexpr int PROJECTION_METHOD_LCC_2SP = 9802;
132 : static_assert(PROJECTION_METHOD_LCC_2SP ==
133 : EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_2SP);
134 : constexpr int PROJECTION_METHOD_OBLIQUE_STEREOGRAPHIC = 9809;
135 : static_assert(PROJECTION_METHOD_OBLIQUE_STEREOGRAPHIC ==
136 : EPSG_CODE_METHOD_OBLIQUE_STEREOGRAPHIC);
137 : constexpr int PROJECTION_METHOD_POLAR_STEREOGRAPHIC = 9810;
138 : static_assert(PROJECTION_METHOD_POLAR_STEREOGRAPHIC ==
139 : EPSG_CODE_METHOD_POLAR_STEREOGRAPHIC_VARIANT_A);
140 : constexpr int PROJECTION_METHOD_KROVAK_OBLIQUE_CONIC_CONFORMAL = 9819;
141 : static_assert(PROJECTION_METHOD_KROVAK_OBLIQUE_CONIC_CONFORMAL ==
142 : EPSG_CODE_METHOD_KROVAK);
143 : constexpr int PROJECTION_METHOD_AMERICAN_POLYCONIC = 9818;
144 : static_assert(PROJECTION_METHOD_AMERICAN_POLYCONIC ==
145 : EPSG_CODE_METHOD_AMERICAN_POLYCONIC);
146 : constexpr int PROJECTION_METHOD_ALBERS_EQUAL_AREA = 9822;
147 : static_assert(PROJECTION_METHOD_ALBERS_EQUAL_AREA ==
148 : EPSG_CODE_METHOD_ALBERS_EQUAL_AREA);
149 : constexpr int PROJECTION_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA = 9820;
150 : static_assert(PROJECTION_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA ==
151 : EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA);
152 :
153 328 : bool S100ReadSRS(const GDALGroup *poRootGroup, OGRSpatialReference &oSRS)
154 : {
155 : // Get SRS
156 328 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
157 984 : auto poHorizontalCRS = poRootGroup->GetAttribute("horizontalCRS");
158 633 : if (poHorizontalCRS &&
159 633 : poHorizontalCRS->GetDataType().GetClass() == GEDTC_NUMERIC)
160 : {
161 : // horizontalCRS is v2.2
162 305 : const int nHorizontalCRS = poHorizontalCRS->ReadAsInt();
163 305 : if (nHorizontalCRS > 0)
164 : {
165 195 : if (oSRS.importFromEPSG(nHorizontalCRS) != OGRERR_NONE)
166 : {
167 0 : oSRS.Clear();
168 : }
169 : }
170 : else
171 : {
172 : auto poNameOfHorizontalCRS =
173 220 : poRootGroup->GetAttribute("nameOfHorizontalCRS");
174 : auto poTypeOfHorizontalCRS =
175 220 : poRootGroup->GetAttribute("typeOfHorizontalCRS");
176 220 : auto poHorizontalCS = poRootGroup->GetAttribute("horizontalCS");
177 : auto poHorizontalDatum =
178 220 : poRootGroup->GetAttribute("horizontalDatum");
179 110 : if (poNameOfHorizontalCRS &&
180 110 : poNameOfHorizontalCRS->GetDataType().GetClass() ==
181 220 : GEDTC_STRING &&
182 220 : poTypeOfHorizontalCRS &&
183 110 : poTypeOfHorizontalCRS->GetDataType().GetClass() ==
184 110 : GEDTC_NUMERIC &&
185 220 : poHorizontalCS &&
186 220 : poHorizontalCS->GetDataType().GetClass() == GEDTC_NUMERIC &&
187 330 : poHorizontalDatum &&
188 110 : poHorizontalDatum->GetDataType().GetClass() == GEDTC_NUMERIC)
189 : {
190 110 : const int nCRSType = poTypeOfHorizontalCRS->ReadAsInt();
191 110 : constexpr int CRS_TYPE_GEOGRAPHIC = 1;
192 110 : constexpr int CRS_TYPE_PROJECTED = 2;
193 :
194 110 : PJ *crs = nullptr;
195 110 : const char *pszCRSName = poNameOfHorizontalCRS->ReadAsString();
196 :
197 110 : const int nHorizontalCS = poHorizontalCS->ReadAsInt();
198 110 : const int nHorizontalDatum = poHorizontalDatum->ReadAsInt();
199 :
200 110 : auto pjCtxt = OSRGetProjTLSContext();
201 110 : if (nHorizontalDatum < 0)
202 : {
203 : auto poNameOfHorizontalDatum =
204 30 : poRootGroup->GetAttribute("nameOfHorizontalDatum");
205 : auto poPrimeMeridian =
206 30 : poRootGroup->GetAttribute("primeMeridian");
207 30 : auto poSpheroid = poRootGroup->GetAttribute("spheroid");
208 15 : if (poNameOfHorizontalDatum &&
209 15 : poNameOfHorizontalDatum->GetDataType().GetClass() ==
210 30 : GEDTC_STRING &&
211 28 : poPrimeMeridian &&
212 13 : poPrimeMeridian->GetDataType().GetClass() ==
213 13 : GEDTC_NUMERIC &&
214 43 : poSpheroid &&
215 13 : poSpheroid->GetDataType().GetClass() == GEDTC_NUMERIC)
216 : {
217 : const char *pszDatumName =
218 13 : poNameOfHorizontalDatum->ReadAsString();
219 13 : const char *pszGeogCRSName =
220 : nCRSType == CRS_TYPE_PROJECTED
221 20 : ? (pszDatumName ? pszDatumName : "unknown")
222 7 : : (pszCRSName ? pszCRSName : "unknown");
223 :
224 13 : const int nPrimeMeridian = poPrimeMeridian->ReadAsInt();
225 13 : const int nSpheroid = poSpheroid->ReadAsInt();
226 13 : PJ *prime_meridian = proj_create_from_database(
227 : pjCtxt, "EPSG", CPLSPrintf("%d", nPrimeMeridian),
228 : PJ_CATEGORY_PRIME_MERIDIAN, false, nullptr);
229 13 : if (!prime_meridian)
230 0 : return false;
231 13 : PJ *spheroid = proj_create_from_database(
232 : pjCtxt, "EPSG", CPLSPrintf("%d", nSpheroid),
233 : PJ_CATEGORY_ELLIPSOID, false, nullptr);
234 13 : if (!spheroid)
235 : {
236 0 : proj_destroy(prime_meridian);
237 0 : return false;
238 : }
239 :
240 13 : double semiMajorMetre = 0;
241 13 : double invFlattening = 0;
242 13 : proj_ellipsoid_get_parameters(pjCtxt, spheroid,
243 : &semiMajorMetre, nullptr,
244 : nullptr, &invFlattening);
245 :
246 13 : double primeMeridianOffset = 0;
247 13 : double primeMeridianUnitConv = 1;
248 13 : const char *primeMeridianUnitName = nullptr;
249 13 : proj_prime_meridian_get_parameters(
250 : pjCtxt, prime_meridian, &primeMeridianOffset,
251 : &primeMeridianUnitConv, &primeMeridianUnitName);
252 :
253 13 : PJ *ellipsoidalCS = proj_create_ellipsoidal_2D_cs(
254 : pjCtxt, PJ_ELLPS2D_LATITUDE_LONGITUDE, nullptr,
255 : 0.0);
256 :
257 13 : crs = proj_create_geographic_crs(
258 : pjCtxt, pszGeogCRSName,
259 : pszDatumName ? pszDatumName : "unknown",
260 : proj_get_name(spheroid), semiMajorMetre,
261 : invFlattening, proj_get_name(prime_meridian),
262 : primeMeridianOffset, primeMeridianUnitName,
263 : primeMeridianUnitConv, ellipsoidalCS);
264 :
265 13 : proj_destroy(ellipsoidalCS);
266 13 : proj_destroy(spheroid);
267 13 : proj_destroy(prime_meridian);
268 : }
269 : else
270 : {
271 2 : CPLError(CE_Warning, CPLE_AppDefined,
272 : "Cannot establish CRS because "
273 : "nameOfHorizontalDatum, primeMeridian and/or "
274 : "spheroid are missing");
275 2 : return false;
276 : }
277 : }
278 : else
279 : {
280 95 : PJ *datum = proj_create_from_database(
281 : pjCtxt, "EPSG", CPLSPrintf("%d", nHorizontalDatum),
282 : PJ_CATEGORY_DATUM, false, nullptr);
283 95 : if (!datum)
284 0 : return false;
285 95 : const char *pszDatumName = proj_get_name(datum);
286 95 : const char *pszGeogCRSName =
287 : nCRSType == CRS_TYPE_PROJECTED
288 105 : ? pszDatumName
289 10 : : (pszCRSName ? pszCRSName : "unknown");
290 :
291 95 : PJ *ellipsoidalCS = proj_create_ellipsoidal_2D_cs(
292 : pjCtxt, PJ_ELLPS2D_LATITUDE_LONGITUDE, nullptr, 0.0);
293 :
294 95 : crs = proj_create_geographic_crs_from_datum(
295 : pjCtxt, pszGeogCRSName, datum, ellipsoidalCS);
296 :
297 95 : proj_destroy(ellipsoidalCS);
298 95 : proj_destroy(datum);
299 : }
300 :
301 108 : if (nCRSType == CRS_TYPE_PROJECTED)
302 : {
303 : auto poProjectionMethod =
304 182 : poRootGroup->GetAttribute("projectionMethod");
305 181 : if (poProjectionMethod &&
306 90 : poProjectionMethod->GetDataType().GetClass() ==
307 91 : GEDTC_NUMERIC)
308 : {
309 : const int nProjectionMethod =
310 90 : poProjectionMethod->ReadAsInt();
311 :
312 : auto poFalseEasting =
313 180 : poRootGroup->GetAttribute("falseEasting");
314 : const double falseEasting =
315 90 : poFalseEasting &&
316 90 : poFalseEasting->GetDataType().GetClass() ==
317 90 : GEDTC_NUMERIC
318 180 : ? poFalseEasting->ReadAsDouble()
319 90 : : 0.0;
320 :
321 : auto poFalseNorthing =
322 180 : poRootGroup->GetAttribute("falseNorthing");
323 : const double falseNorthing =
324 90 : poFalseNorthing &&
325 90 : poFalseNorthing->GetDataType().GetClass() ==
326 90 : GEDTC_NUMERIC
327 180 : ? poFalseNorthing->ReadAsDouble()
328 90 : : 0.0;
329 :
330 299 : const auto getParameter = [poRootGroup](int nParam)
331 : {
332 : auto poParam = poRootGroup->GetAttribute(
333 897 : "projectionParameter" + std::to_string(nParam));
334 299 : if (poParam && poParam->GetDataType().GetClass() ==
335 299 : GEDTC_NUMERIC)
336 299 : return poParam->ReadAsDouble();
337 0 : CPLError(
338 : CE_Warning, CPLE_AppDefined, "%s",
339 0 : std::string(
340 : "Missing attribute projectionParameter")
341 0 : .append(std::to_string(nParam))
342 : .c_str());
343 0 : return 0.0;
344 90 : };
345 :
346 90 : PJ *conv = nullptr;
347 90 : switch (nProjectionMethod)
348 : {
349 13 : case PROJECTION_METHOD_MERCATOR:
350 : {
351 : conv =
352 13 : proj_create_conversion_mercator_variant_b(
353 : pjCtxt,
354 : /* latitude_first_parallel = */
355 : getParameter(1),
356 : /* center_long = */ getParameter(2),
357 : falseEasting, falseNorthing, nullptr,
358 : 0.0, nullptr, 0.0);
359 13 : break;
360 : }
361 :
362 7 : case PROJECTION_METHOD_TRANSVERSE_MERCATOR:
363 : {
364 : conv =
365 7 : proj_create_conversion_transverse_mercator(
366 : pjCtxt,
367 : /* center_lat = */ getParameter(1),
368 : /* center_long = */ getParameter(2),
369 : /* scale = */ getParameter(3),
370 : falseEasting, falseNorthing, nullptr,
371 : 0.0, nullptr, 0.0);
372 7 : break;
373 : }
374 :
375 7 : case PROJECTION_METHOD_OBLIQUE_MERCATOR:
376 : {
377 : conv =
378 7 : proj_create_conversion_hotine_oblique_mercator_variant_b(
379 : pjCtxt,
380 : /* latitude_projection_centre = */
381 : getParameter(1),
382 : /* longitude_projection_centre = */
383 : getParameter(2),
384 : /* azimuth_initial_line = */
385 : getParameter(3),
386 : /* angle_from_rectified_to_skrew_grid = */
387 : getParameter(4),
388 : /* scale = */ getParameter(5),
389 : falseEasting, falseNorthing, nullptr,
390 : 0.0, nullptr, 0.0);
391 7 : break;
392 : }
393 :
394 7 : case PROJECTION_METHOD_HOTINE_OBLIQUE_MERCATOR:
395 : {
396 : conv =
397 7 : proj_create_conversion_hotine_oblique_mercator_variant_a(
398 : pjCtxt,
399 : /* latitude_projection_centre = */
400 : getParameter(1),
401 : /* longitude_projection_centre = */
402 : getParameter(2),
403 : /* azimuth_initial_line = */
404 : getParameter(3),
405 : /* angle_from_rectified_to_skrew_grid = */
406 : getParameter(4),
407 : /* scale = */ getParameter(5),
408 : falseEasting, falseNorthing, nullptr,
409 : 0.0, nullptr, 0.0);
410 7 : break;
411 : }
412 :
413 7 : case PROJECTION_METHOD_LCC_1SP:
414 : {
415 : conv =
416 7 : proj_create_conversion_lambert_conic_conformal_1sp(
417 : pjCtxt,
418 : /* center_lat = */ getParameter(1),
419 : /* center_long = */ getParameter(2),
420 : /* scale = */ getParameter(3),
421 : falseEasting, falseNorthing, nullptr,
422 : 0.0, nullptr, 0.0);
423 7 : break;
424 : }
425 :
426 7 : case PROJECTION_METHOD_LCC_2SP:
427 : {
428 : conv =
429 7 : proj_create_conversion_lambert_conic_conformal_2sp(
430 : pjCtxt,
431 : /* latitude_false_origin = */
432 : getParameter(1),
433 : /* longitude_false_origin = */
434 : getParameter(2),
435 : /* latitude_first_parallel = */
436 : getParameter(3),
437 : /* latitude_second_parallel = */
438 : getParameter(4), falseEasting,
439 : falseNorthing, nullptr, 0.0, nullptr,
440 : 0.0);
441 7 : break;
442 : }
443 :
444 7 : case PROJECTION_METHOD_OBLIQUE_STEREOGRAPHIC:
445 : {
446 : conv =
447 7 : proj_create_conversion_oblique_stereographic(
448 : pjCtxt,
449 : /* center_lat = */ getParameter(1),
450 : /* center_long = */ getParameter(2),
451 : /* scale = */ getParameter(3),
452 : falseEasting, falseNorthing, nullptr,
453 : 0.0, nullptr, 0.0);
454 7 : break;
455 : }
456 :
457 7 : case PROJECTION_METHOD_POLAR_STEREOGRAPHIC:
458 : {
459 : conv =
460 7 : proj_create_conversion_polar_stereographic_variant_a(
461 : pjCtxt,
462 : /* center_lat = */ getParameter(1),
463 : /* center_long = */ getParameter(2),
464 : /* scale = */ getParameter(3),
465 : falseEasting, falseNorthing, nullptr,
466 : 0.0, nullptr, 0.0);
467 7 : break;
468 : }
469 :
470 7 : case PROJECTION_METHOD_KROVAK_OBLIQUE_CONIC_CONFORMAL:
471 : {
472 7 : conv = proj_create_conversion_krovak(
473 : pjCtxt,
474 : /* latitude_projection_centre = */
475 : getParameter(1),
476 : /* longitude_of_origin = */ getParameter(2),
477 : /* colatitude_cone_axis = */
478 : getParameter(3),
479 : /* latitude_pseudo_standard_parallel = */
480 : getParameter(4),
481 : /* scale_factor_pseudo_standard_parallel = */
482 : getParameter(5), falseEasting,
483 : falseNorthing, nullptr, 0.0, nullptr, 0.0);
484 7 : break;
485 : }
486 :
487 7 : case PROJECTION_METHOD_AMERICAN_POLYCONIC:
488 : {
489 : conv =
490 7 : proj_create_conversion_american_polyconic(
491 : pjCtxt,
492 : /* center_lat = */ getParameter(1),
493 : /* center_long = */ getParameter(2),
494 : falseEasting, falseNorthing, nullptr,
495 : 0.0, nullptr, 0.0);
496 7 : break;
497 : }
498 :
499 7 : case PROJECTION_METHOD_ALBERS_EQUAL_AREA:
500 : {
501 7 : conv = proj_create_conversion_albers_equal_area(
502 : pjCtxt,
503 : /* latitude_false_origin = */
504 : getParameter(1),
505 : /* longitude_false_origin = */
506 : getParameter(2),
507 : /* latitude_first_parallel = */
508 : getParameter(3),
509 : /* latitude_second_parallel = */
510 : getParameter(4), falseEasting,
511 : falseNorthing, nullptr, 0.0, nullptr, 0.0);
512 7 : break;
513 : }
514 :
515 7 : case PROJECTION_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA:
516 : {
517 : conv =
518 7 : proj_create_conversion_lambert_azimuthal_equal_area(
519 : pjCtxt,
520 : /* latitude_nat_origin = */
521 : getParameter(1),
522 : /* longitude_nat_origin = */
523 : getParameter(2), falseEasting,
524 : falseNorthing, nullptr, 0.0, nullptr,
525 : 0.0);
526 7 : break;
527 : }
528 :
529 0 : default:
530 : {
531 0 : CPLError(CE_Warning, CPLE_AppDefined,
532 : "Cannot establish CRS because "
533 : "projectionMethod = %d is unhandled",
534 : nProjectionMethod);
535 0 : proj_destroy(crs);
536 0 : return false;
537 : }
538 : }
539 :
540 90 : constexpr int HORIZONTAL_CS_EASTING_NORTHING_METRE =
541 : 4400;
542 : // constexpr int HORIZONTAL_CS_NORTHING_EASTING_METRE = 4500;
543 180 : PJ *cartesianCS = proj_create_cartesian_2D_cs(
544 : pjCtxt,
545 : nHorizontalCS ==
546 : HORIZONTAL_CS_EASTING_NORTHING_METRE
547 90 : ? PJ_CART2D_EASTING_NORTHING
548 : : PJ_CART2D_NORTHING_EASTING,
549 : nullptr, 0.0);
550 :
551 90 : PJ *projCRS = proj_create_projected_crs(
552 : pjCtxt, pszCRSName ? pszCRSName : "unknown", crs,
553 : conv, cartesianCS);
554 90 : proj_destroy(crs);
555 90 : crs = projCRS;
556 90 : proj_destroy(conv);
557 90 : proj_destroy(cartesianCS);
558 : }
559 : else
560 : {
561 1 : CPLError(CE_Warning, CPLE_AppDefined,
562 : "Cannot establish CRS because "
563 : "projectionMethod is missing");
564 1 : proj_destroy(crs);
565 1 : return false;
566 : }
567 : }
568 17 : else if (nCRSType != CRS_TYPE_GEOGRAPHIC)
569 : {
570 0 : CPLError(CE_Warning, CPLE_AppDefined,
571 : "Cannot establish CRS because of invalid value "
572 : "for typeOfHorizontalCRS: %d",
573 : nCRSType);
574 0 : proj_destroy(crs);
575 0 : return false;
576 : }
577 :
578 : const char *pszPROJJSON =
579 107 : proj_as_projjson(pjCtxt, crs, nullptr);
580 107 : oSRS.SetFromUserInput(pszPROJJSON);
581 :
582 107 : proj_destroy(crs);
583 : }
584 : else
585 : {
586 0 : CPLError(CE_Warning, CPLE_AppDefined,
587 : "Cannot establish CRS because nameOfHorizontalCRS, "
588 : "typeOfHorizontalCRS, horizontalCS and/or "
589 : "horizontalDatum are missing");
590 0 : return false;
591 : }
592 : }
593 : }
594 : else
595 : {
596 : auto poHorizontalDatumReference =
597 69 : poRootGroup->GetAttribute("horizontalDatumReference");
598 : auto poHorizontalDatumValue =
599 69 : poRootGroup->GetAttribute("horizontalDatumValue");
600 23 : if (poHorizontalDatumReference && poHorizontalDatumValue)
601 : {
602 : const char *pszAuthName =
603 23 : poHorizontalDatumReference->ReadAsString();
604 23 : const char *pszAuthCode = poHorizontalDatumValue->ReadAsString();
605 23 : if (pszAuthName && pszAuthCode)
606 : {
607 46 : if (oSRS.SetFromUserInput(
608 46 : (std::string(pszAuthName) + ':' + pszAuthCode).c_str(),
609 : OGRSpatialReference::
610 23 : SET_FROM_USER_INPUT_LIMITATIONS_get()) !=
611 : OGRERR_NONE)
612 : {
613 0 : oSRS.Clear();
614 : }
615 : }
616 : }
617 : }
618 325 : return !oSRS.IsEmpty();
619 : }
620 :
621 : /************************************************************************/
622 : /* S100GetNumPointsLongitudinalLatitudinal() */
623 : /************************************************************************/
624 :
625 454 : bool S100GetNumPointsLongitudinalLatitudinal(const GDALGroup *poGroup,
626 : int &nNumPointsLongitudinal,
627 : int &nNumPointsLatitudinal)
628 : {
629 1362 : auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
630 1362 : auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
631 : auto poNumPointsLongitudinal =
632 1362 : poGroup->GetAttribute("numPointsLongitudinal");
633 1362 : auto poNumPointsLatitudinal = poGroup->GetAttribute("numPointsLatitudinal");
634 454 : if (poSpacingX &&
635 1362 : poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
636 908 : poSpacingY &&
637 908 : poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64 &&
638 454 : poNumPointsLongitudinal &&
639 454 : GDALDataTypeIsInteger(
640 908 : poNumPointsLongitudinal->GetDataType().GetNumericDataType()) &&
641 1362 : poNumPointsLatitudinal &&
642 454 : GDALDataTypeIsInteger(
643 454 : poNumPointsLatitudinal->GetDataType().GetNumericDataType()))
644 : {
645 454 : nNumPointsLongitudinal = poNumPointsLongitudinal->ReadAsInt();
646 454 : nNumPointsLatitudinal = poNumPointsLatitudinal->ReadAsInt();
647 :
648 : // Those are optional, but use them when available, to detect
649 : // potential inconsistency
650 1362 : auto poEastBoundLongitude = poGroup->GetAttribute("eastBoundLongitude");
651 1362 : auto poWestBoundLongitude = poGroup->GetAttribute("westBoundLongitude");
652 : auto poSouthBoundLongitude =
653 1362 : poGroup->GetAttribute("southBoundLatitude");
654 908 : auto poNorthBoundLatitude = poGroup->GetAttribute("northBoundLatitude");
655 807 : if (poEastBoundLongitude &&
656 353 : GDALDataTypeIsFloating(
657 1160 : poEastBoundLongitude->GetDataType().GetNumericDataType()) &&
658 353 : poWestBoundLongitude &&
659 353 : GDALDataTypeIsFloating(
660 706 : poWestBoundLongitude->GetDataType().GetNumericDataType()) &&
661 353 : poSouthBoundLongitude &&
662 353 : GDALDataTypeIsFloating(
663 706 : poSouthBoundLongitude->GetDataType().GetNumericDataType()) &&
664 1160 : poNorthBoundLatitude &&
665 353 : GDALDataTypeIsFloating(
666 353 : poNorthBoundLatitude->GetDataType().GetNumericDataType()))
667 : {
668 353 : const double dfSpacingX = poSpacingX->ReadAsDouble();
669 353 : const double dfSpacingY = poSpacingY->ReadAsDouble();
670 :
671 353 : const double dfEast = poEastBoundLongitude->ReadAsDouble();
672 353 : const double dfWest = poWestBoundLongitude->ReadAsDouble();
673 353 : const double dfSouth = poSouthBoundLongitude->ReadAsDouble();
674 353 : const double dfNorth = poNorthBoundLatitude->ReadAsDouble();
675 353 : if (std::fabs((dfWest + nNumPointsLongitudinal * dfSpacingX) -
676 353 : dfEast) < 5 * dfSpacingX &&
677 353 : std::fabs((dfSouth + nNumPointsLatitudinal * dfSpacingY) -
678 353 : dfNorth) < 5 * dfSpacingY)
679 : {
680 : // We need up to 5 spacings for product
681 : // S-111 Trial Data Set Release 1.1/111UK_20210401T000000Z_SolentAndAppr_dcf2.h5
682 : }
683 : else
684 : {
685 0 : CPLError(
686 : CE_Warning, CPLE_AppDefined,
687 : "Caution: "
688 : "eastBoundLongitude/westBoundLongitude/southBoundLatitude/"
689 : "northBoundLatitude/gridSpacingLongitudinal/"
690 : "gridSpacingLatitudinal/numPointsLongitudinal/"
691 : "numPointsLatitudinal do not seem to be consistent");
692 0 : CPLDebug("S100", "Computed east = %f. Actual = %f",
693 0 : dfWest + nNumPointsLongitudinal * dfSpacingX, dfEast);
694 0 : CPLDebug("S100", "Computed north = %f. Actual = %f",
695 0 : dfSouth + nNumPointsLatitudinal * dfSpacingY, dfNorth);
696 : }
697 : }
698 :
699 454 : return true;
700 : }
701 0 : return false;
702 : }
703 :
704 : /************************************************************************/
705 : /* S100GetGeoTransform() */
706 : /************************************************************************/
707 :
708 197 : bool S100GetGeoTransform(const GDALGroup *poGroup, GDALGeoTransform >,
709 : bool bNorthUp)
710 : {
711 591 : auto poOriginX = poGroup->GetAttribute("gridOriginLongitude");
712 591 : auto poOriginY = poGroup->GetAttribute("gridOriginLatitude");
713 591 : auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
714 591 : auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
715 197 : if (poOriginX &&
716 591 : poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 &&
717 394 : poOriginY &&
718 394 : poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 &&
719 394 : poSpacingX &&
720 394 : poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
721 591 : poSpacingY &&
722 197 : poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64)
723 : {
724 197 : int nNumPointsLongitudinal = 0;
725 197 : int nNumPointsLatitudinal = 0;
726 197 : if (!S100GetNumPointsLongitudinalLatitudinal(
727 : poGroup, nNumPointsLongitudinal, nNumPointsLatitudinal))
728 0 : return false;
729 :
730 197 : const double dfSpacingX = poSpacingX->ReadAsDouble();
731 197 : const double dfSpacingY = poSpacingY->ReadAsDouble();
732 :
733 197 : gt.xorig = poOriginX->ReadAsDouble();
734 197 : gt.yorig = poOriginY->ReadAsDouble() +
735 197 : (bNorthUp ? dfSpacingY * (nNumPointsLatitudinal - 1) : 0);
736 197 : gt.xscale = dfSpacingX;
737 197 : gt.yscale = bNorthUp ? -dfSpacingY : dfSpacingY;
738 :
739 : // From pixel-center convention to pixel-corner convention
740 197 : gt.xorig -= gt.xscale / 2;
741 197 : gt.yorig -= gt.yscale / 2;
742 :
743 197 : return true;
744 : }
745 0 : return false;
746 : }
747 :
748 : /************************************************************************/
749 : /* S100GetDimensions() */
750 : /************************************************************************/
751 :
752 104 : bool S100GetDimensions(
753 : const GDALGroup *poGroup,
754 : std::vector<std::shared_ptr<GDALDimension>> &apoDims,
755 : std::vector<std::shared_ptr<GDALMDArray>> &apoIndexingVars)
756 : {
757 312 : auto poOriginX = poGroup->GetAttribute("gridOriginLongitude");
758 312 : auto poOriginY = poGroup->GetAttribute("gridOriginLatitude");
759 312 : auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
760 312 : auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
761 104 : if (poOriginX &&
762 312 : poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 &&
763 208 : poOriginY &&
764 208 : poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 &&
765 208 : poSpacingX &&
766 208 : poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
767 312 : poSpacingY &&
768 104 : poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64)
769 : {
770 104 : int nNumPointsLongitudinal = 0;
771 104 : int nNumPointsLatitudinal = 0;
772 104 : if (!S100GetNumPointsLongitudinalLatitudinal(
773 : poGroup, nNumPointsLongitudinal, nNumPointsLatitudinal))
774 0 : return false;
775 :
776 : {
777 : auto poDim = std::make_shared<GDALDimensionWeakIndexingVar>(
778 208 : std::string(), "Y", GDAL_DIM_TYPE_HORIZONTAL_Y, std::string(),
779 208 : nNumPointsLatitudinal);
780 : auto poIndexingVar = GDALMDArrayRegularlySpaced::Create(
781 312 : std::string(), poDim->GetName(), poDim,
782 416 : poOriginY->ReadAsDouble(), poSpacingY->ReadAsDouble(), 0);
783 104 : poDim->SetIndexingVariable(poIndexingVar);
784 104 : apoDims.emplace_back(poDim);
785 104 : apoIndexingVars.emplace_back(poIndexingVar);
786 : }
787 :
788 : {
789 : auto poDim = std::make_shared<GDALDimensionWeakIndexingVar>(
790 208 : std::string(), "X", GDAL_DIM_TYPE_HORIZONTAL_X, std::string(),
791 208 : nNumPointsLongitudinal);
792 : auto poIndexingVar = GDALMDArrayRegularlySpaced::Create(
793 312 : std::string(), poDim->GetName(), poDim,
794 416 : poOriginX->ReadAsDouble(), poSpacingX->ReadAsDouble(), 0);
795 104 : poDim->SetIndexingVariable(poIndexingVar);
796 104 : apoDims.emplace_back(poDim);
797 104 : apoIndexingVars.emplace_back(poIndexingVar);
798 : }
799 :
800 104 : return true;
801 : }
802 0 : return false;
803 : }
804 :
805 : /************************************************************************/
806 : /* SetMetadataForDataDynamicity() */
807 : /************************************************************************/
808 :
809 123 : void S100BaseDataset::SetMetadataForDataDynamicity(const GDALAttribute *poAttr)
810 : {
811 123 : if (poAttr && poAttr->GetDataType().GetClass() == GEDTC_NUMERIC)
812 : {
813 123 : const int nVal = poAttr->ReadAsInt();
814 : const std::map<int, std::pair<const char *, const char *>> map = {
815 : {1,
816 : {"observation", "Values from in-situ sensor(s); may be quality "
817 : "controlled and stored after collection"}},
818 : {2,
819 : {"astronomicalPrediction",
820 : "Values computed using harmonic analysis or other proven method "
821 : "of tidal analysis"}},
822 : {3,
823 : {"analysisOrHybrid",
824 : "Values calculated by statistical or other indirect methods, or "
825 : "a combination of methods"}},
826 : {3,
827 : {"hydrodynamicHindcast",
828 : "Values calculated from a two- or three-dimensional dynamic "
829 : "simulation of past conditions using only observed data for "
830 : "boundary forcing, via statistical method or combination"}},
831 : {5,
832 : {"hydrodynamicForecast",
833 : "Values calculated from a two- or three-dimensional dynamic "
834 : "simulation of future conditions using predicted data for "
835 : "boundary forcing, via statistical method or combination"}},
836 : {6,
837 : {"observedMinusPredicted",
838 : "Values computed as observed minus predicted values"}},
839 : {7,
840 : {"observedMinusAnalysis",
841 : "Values computed as observed minus analysis values"}},
842 : {8,
843 : {"observedMinusHindcast",
844 : "Values computed as observed minus hindcast values"}},
845 : {9,
846 : {"observedMinusForecast",
847 : "Values computed as observed minus forecast values"}},
848 : {10,
849 : {"forecastMinusPredicted",
850 : "Values computed as forecast minus predicted values"}},
851 246 : };
852 123 : const auto oIter = map.find(nVal);
853 123 : if (oIter != map.end())
854 : {
855 246 : GDALDataset::SetMetadataItem("DATA_DYNAMICITY_NAME",
856 123 : oIter->second.first);
857 246 : GDALDataset::SetMetadataItem("DATA_DYNAMICITY_DEFINITION",
858 123 : oIter->second.second);
859 : }
860 : }
861 123 : }
862 :
863 : /************************************************************************/
864 : /* SetMetadataForCommonPointRule() */
865 : /************************************************************************/
866 :
867 139 : void S100BaseDataset::SetMetadataForCommonPointRule(const GDALAttribute *poAttr)
868 : {
869 139 : if (poAttr && poAttr->GetDataType().GetClass() == GEDTC_NUMERIC)
870 : {
871 139 : const int nVal = poAttr->ReadAsInt();
872 : const std::map<int, std::pair<const char *, const char *>> map = {
873 : {1, {"average", "return the mean of the attribute values"}},
874 : {2, {"low", "use the least of the attribute values"}},
875 : {3, {"high", "use the greatest of the attribute values"}},
876 : {4,
877 : {"all", "return all the attribute values that can be determined "
878 : "for the position"}},
879 278 : };
880 139 : const auto oIter = map.find(nVal);
881 139 : if (oIter != map.end())
882 : {
883 278 : GDALDataset::SetMetadataItem("COMMON_POINT_RULE_NAME",
884 139 : oIter->second.first);
885 278 : GDALDataset::SetMetadataItem("COMMON_POINT_RULE_DEFINITION",
886 139 : oIter->second.second);
887 : }
888 : }
889 139 : }
890 :
891 : /************************************************************************/
892 : /* SetMetadataForInterpolationType() */
893 : /************************************************************************/
894 :
895 139 : void S100BaseDataset::SetMetadataForInterpolationType(
896 : const GDALAttribute *poAttr)
897 : {
898 139 : if (poAttr && poAttr->GetDataType().GetClass() == GEDTC_NUMERIC)
899 : {
900 139 : const int nVal = poAttr->ReadAsInt();
901 : const std::map<int, std::pair<const char *, const char *>> map = {
902 : {1,
903 : {"nearestneighbor",
904 : "Assign the feature attribute value associated with the nearest "
905 : "domain object in the domain of the coverage"}},
906 : {5,
907 : {"bilinear", "Assign a value computed by using a bilinear "
908 : "function of position within the grid cell"}},
909 : {6,
910 : {"biquadratic", "Assign a value computed by using a biquadratic "
911 : "function of position within the grid cell"}},
912 : {7,
913 : {"bicubic", "Assign a value computed by using a bicubic function "
914 : "of position within the grid cell"}},
915 : {8,
916 : {"lostarea", "Assign a value computed by using the lost area "
917 : "method described in ISO 19123"}},
918 : {9,
919 : {"barycentric", "Assign a value computed by using the barycentric "
920 : "method described in ISO 19123"}},
921 : {10,
922 : {"discrete", "No interpolation method applies to the coverage"}},
923 278 : };
924 139 : const auto oIter = map.find(nVal);
925 139 : if (oIter != map.end())
926 : {
927 278 : GDALDataset::SetMetadataItem("INTERPOLATION_TYPE_NAME",
928 139 : oIter->second.first);
929 278 : GDALDataset::SetMetadataItem("INTERPOLATION_TYPE_DEFINITION",
930 139 : oIter->second.second);
931 : }
932 : }
933 139 : }
934 :
935 : /************************************************************************/
936 : /* gasVerticalDatums */
937 : /************************************************************************/
938 :
939 : // https://iho.int/uploads/user/pubs/standards/s-100/S-100_5.2.0_Final_Clean.pdf
940 : // Table 10c-25 - Vertical and sounding datum, page 53
941 : static const struct
942 : {
943 : int nCode;
944 : const char *pszName;
945 : const char *pszAbbrev;
946 : const char *pszDefinition;
947 : } gasVerticalDatums[] = {
948 : {1, "meanLowWaterSprings", "MLWS",
949 : "The average height of the low waters of spring tides. This level is used "
950 : "as a tidal datum in some areas. Also called spring low water."},
951 : {2, "meanLowerLowWaterSprings", nullptr,
952 : "The average height of lower low water springs at a place"},
953 : {3, "meanSeaLevel", "MSL",
954 : "The average height of the surface of the sea at a tide station for all "
955 : "stages of the tide over a 19-year period, usually determined from hourly "
956 : "height readings measured from a fixed predetermined reference level."},
957 : {4, "lowestLowWater", nullptr,
958 : "An arbitrary level conforming to the lowest tide observed at a place, or "
959 : "some what lower."},
960 : {5, "meanLowWater", "MLW",
961 : "The average height of all low waters at a place over a 19-year period."},
962 : {6, "lowestLowWaterSprings", nullptr,
963 : "An arbitrary level conforming to the lowest water level observed at a "
964 : "place at spring tides during a period of time shorter than 19 years."},
965 : {7, "approximateMeanLowWaterSprings", nullptr,
966 : "An arbitrary level, usually within 0.3m from that of Mean Low Water "
967 : "Springs (MLWS)."},
968 : {8, "indianSpringLowWater", "ISLW",
969 : "An arbitrary tidal datum approximating the level of the mean of the "
970 : "lower low water at spring tides. It was first used in waters surrounding "
971 : "India."},
972 : {9, "lowWaterSprings", nullptr,
973 : "An arbitrary level, approximating that of mean low water springs "
974 : "(MLWS)."},
975 : {10, "approximateLowestAstronomicalTide", nullptr,
976 : "An arbitrary level, usually within 0.3m from that of Lowest Astronomical "
977 : "Tide (LAT)."},
978 : {11, "nearlyLowestLowWater", nullptr,
979 : "An arbitrary level approximating the lowest water level observed at a "
980 : "place, usually equivalent to the Indian Spring Low Water (ISLW)."},
981 : {12, "meanLowerLowWater", "MLLW",
982 : "The average height of the lower low waters at a place over a 19-year "
983 : "period."},
984 : {13, "lowWater", "LW",
985 : "The lowest level reached at a place by the water surface in one "
986 : "oscillation. Also called low tide."},
987 : {14, "approximateMeanLowWater", nullptr,
988 : "An arbitrary level, usually within 0.3m from that of Mean Low Water "
989 : "(MLW)."},
990 : {15, "approximateMeanLowerLowWater", nullptr,
991 : "An arbitrary level, usually within 0.3m from that of Mean Lower Low "
992 : "Water (MLLW)."},
993 : {16, "meanHighWater", "MHW",
994 : "The average height of all high waters at a place over a 19-year period."},
995 : {17, "meanHighWaterSprings", "MHWS",
996 : "The average height of the high waters of spring tides. Also called "
997 : "spring high water."},
998 : {18, "highWater", "HW",
999 : "The highest level reached at a place by the water surface in one "
1000 : "oscillation."},
1001 : {19, "approximateMeanSeaLevel", nullptr,
1002 : "An arbitrary level, usually within 0.3m from that of Mean Sea Level "
1003 : "(MSL)."},
1004 : {20, "highWaterSprings", nullptr,
1005 : "An arbitrary level, approximating that of mean high water springs "
1006 : "(MHWS)."},
1007 : {21, "meanHigherHighWater", "MHHW",
1008 : "The average height of higher high waters at a place over a 19-year "
1009 : "period."},
1010 : {22, "equinoctialSpringLowWater", nullptr,
1011 : "The level of low water springs near the time of an equinox."},
1012 : {23, "lowestAstronomicalTide", "LAT",
1013 : "The lowest tide level which can be predicted to occur under average "
1014 : "meteorological conditions and under any combination of astronomical "
1015 : "conditions."},
1016 : {24, "localDatum", nullptr,
1017 : "An arbitrary datum defined by a local harbour authority, from which "
1018 : "levels and tidal heights are measured by this authority."},
1019 : {25, "internationalGreatLakesDatum1985", nullptr,
1020 : "A vertical reference system with its zero based on the mean water level "
1021 : "at Rimouski/Pointe-au-Pere, Quebec, over the period 1970 to 1988."},
1022 : {26, "meanWaterLevel", nullptr,
1023 : "The average of all hourly water levels over the available period of "
1024 : "record."},
1025 : {27, "lowerLowWaterLargeTide", nullptr,
1026 : "The average of the lowest low waters, one from each of 19 years of "
1027 : "observations."},
1028 : {28, "higherHighWaterLargeTide", nullptr,
1029 : "The average of the highest high waters, one from each of 19 years of "
1030 : "observations."},
1031 : {29, "nearlyHighestHighWater", nullptr,
1032 : "An arbitrary level approximating the highest water level observed at a "
1033 : "place, usually equivalent to the high water springs."},
1034 : {30, "highestAstronomicalTide", "HAT",
1035 : "The highest tidal level which can be predicted to occur under average "
1036 : "meteorological conditions and under any combination of astronomical "
1037 : "conditions."},
1038 : {44, "balticSeaChartDatum2000", nullptr,
1039 : "The datum refers to each Baltic country's realization of the European "
1040 : "Vertical Reference System (EVRS) with land-uplift epoch 2000, which is "
1041 : "connected to the Normaal Amsterdams Peil (NAP)"},
1042 : {46, "internationalGreatLakesDatum2020", nullptr,
1043 : "The 2020 update to the International Great Lakes Datum, the official "
1044 : "reference system used to measure water level heights in the Great Lakes, "
1045 : "connecting channels, and the St. Lawrence River system."},
1046 : {47, "seaFloor", nullptr,
1047 : "The bottom of the ocean and seas where there is a generally smooth "
1048 : "gentle gradient. Also referred to as sea bed (sometimes seabed or "
1049 : "sea-bed), and sea bottom."},
1050 : {48, "seaSurface", nullptr,
1051 : "A two-dimensional (in the horizontal plane) field representing the "
1052 : "air-sea interface, with high-frequency fluctuations such as wind waves "
1053 : "and swell, but not astronomical tides, filtered out"},
1054 : {49, "hydrographicZero", nullptr,
1055 : "A vertical reference near the lowest astronomical tide (LAT) below which "
1056 : "the sea level falls only very exceptionally"},
1057 : };
1058 :
1059 : /************************************************************************/
1060 : /* S100GetVerticalDatumCodeFromNameOrAbbrev() */
1061 : /************************************************************************/
1062 :
1063 88 : int S100GetVerticalDatumCodeFromNameOrAbbrev(const char *pszStr)
1064 : {
1065 88 : const int nCode = atoi(pszStr);
1066 1145 : for (const auto &sEntry : gasVerticalDatums)
1067 : {
1068 1141 : if (sEntry.nCode == nCode || EQUAL(pszStr, sEntry.pszName) ||
1069 1139 : (sEntry.pszAbbrev && EQUAL(pszStr, sEntry.pszAbbrev)))
1070 : {
1071 84 : return sEntry.nCode;
1072 : }
1073 : }
1074 4 : return -1;
1075 : }
1076 :
1077 : /************************************************************************/
1078 : /* S100ReadVerticalDatum() */
1079 : /************************************************************************/
1080 :
1081 454 : void S100ReadVerticalDatum(GDALMajorObject *poMO, const GDALGroup *poGroup)
1082 : {
1083 :
1084 454 : int nVerticalDatumReference = 1;
1085 : auto poVerticalDatumReference =
1086 1362 : poGroup->GetAttribute("verticalDatumReference");
1087 626 : if (poVerticalDatumReference &&
1088 626 : poVerticalDatumReference->GetDataType().GetClass() == GEDTC_NUMERIC)
1089 : {
1090 172 : nVerticalDatumReference = poVerticalDatumReference->ReadAsInt();
1091 : }
1092 :
1093 1362 : auto poVerticalDatum = poGroup->GetAttribute("verticalDatum");
1094 637 : if (poVerticalDatum &&
1095 637 : poVerticalDatum->GetDataType().GetClass() == GEDTC_NUMERIC)
1096 : {
1097 183 : poMO->GDALMajorObject::SetMetadataItem(S100_VERTICAL_DATUM_ABBREV,
1098 : nullptr);
1099 183 : poMO->GDALMajorObject::SetMetadataItem(S100_VERTICAL_DATUM_EPSG_CODE,
1100 : nullptr);
1101 183 : poMO->GDALMajorObject::SetMetadataItem(S100_VERTICAL_DATUM_NAME,
1102 : nullptr);
1103 183 : poMO->GDALMajorObject::SetMetadataItem("verticalDatum", nullptr);
1104 183 : if (nVerticalDatumReference == 1)
1105 : {
1106 177 : bool bFound = false;
1107 177 : const auto nVal = poVerticalDatum->ReadAsInt();
1108 2116 : for (const auto &sVerticalDatum : gasVerticalDatums)
1109 : {
1110 2116 : if (sVerticalDatum.nCode == nVal)
1111 : {
1112 177 : bFound = true;
1113 177 : poMO->GDALMajorObject::SetMetadataItem(
1114 177 : S100_VERTICAL_DATUM_NAME, sVerticalDatum.pszName);
1115 177 : if (sVerticalDatum.pszAbbrev)
1116 : {
1117 177 : poMO->GDALMajorObject::SetMetadataItem(
1118 : S100_VERTICAL_DATUM_ABBREV,
1119 177 : sVerticalDatum.pszAbbrev);
1120 : }
1121 177 : poMO->GDALMajorObject::SetMetadataItem(
1122 : "VERTICAL_DATUM_DEFINITION",
1123 177 : sVerticalDatum.pszDefinition);
1124 177 : break;
1125 : }
1126 : }
1127 177 : if (!bFound)
1128 : {
1129 0 : poMO->GDALMajorObject::SetMetadataItem("verticalDatum",
1130 : CPLSPrintf("%d", nVal));
1131 : }
1132 : }
1133 6 : else if (nVerticalDatumReference == 2)
1134 : {
1135 6 : const auto nVal = poVerticalDatum->ReadAsInt();
1136 6 : PJ *datum = proj_create_from_database(
1137 : OSRGetProjTLSContext(), "EPSG", CPLSPrintf("%d", nVal),
1138 : PJ_CATEGORY_DATUM, false, nullptr);
1139 6 : poMO->GDALMajorObject::SetMetadataItem(
1140 : S100_VERTICAL_DATUM_EPSG_CODE, CPLSPrintf("%d", nVal));
1141 6 : if (datum)
1142 : {
1143 6 : poMO->GDALMajorObject::SetMetadataItem(S100_VERTICAL_DATUM_NAME,
1144 : proj_get_name(datum));
1145 6 : proj_destroy(datum);
1146 : }
1147 : }
1148 : else
1149 : {
1150 0 : CPLDebug("S100", "Unhandled verticalDatumReference = %d",
1151 : nVerticalDatumReference);
1152 : }
1153 : }
1154 454 : }
1155 :
1156 : /************************************************************************/
1157 : /* S100ReadMetadata() */
1158 : /************************************************************************/
1159 :
1160 224 : std::string S100ReadMetadata(GDALDataset *poDS, const std::string &osFilename,
1161 : const GDALGroup *poRootGroup)
1162 : {
1163 224 : std::string osMetadataFile;
1164 3527 : for (const auto &poAttr : poRootGroup->GetAttributes())
1165 : {
1166 3303 : const auto &osName = poAttr->GetName();
1167 3303 : if (osName == "metadata")
1168 : {
1169 30 : const char *pszVal = poAttr->ReadAsString();
1170 30 : if (pszVal && pszVal[0])
1171 : {
1172 30 : if (CPLHasPathTraversal(pszVal))
1173 : {
1174 0 : CPLError(CE_Warning, CPLE_AppDefined,
1175 : "Path traversal detected in %s", pszVal);
1176 0 : continue;
1177 : }
1178 60 : osMetadataFile = CPLFormFilenameSafe(
1179 60 : CPLGetPathSafe(osFilename.c_str()).c_str(), pszVal,
1180 30 : nullptr);
1181 : VSIStatBufL sStat;
1182 30 : if (VSIStatL(osMetadataFile.c_str(), &sStat) != 0)
1183 : {
1184 : // Test products from https://data.admiralty.co.uk/portal/apps/sites/#/marine-data-portal/pages/s-100
1185 : // advertise a metadata filename starting with "MD_", per the spec,
1186 : // but the actual filename does not start with "MD_"...
1187 6 : if (STARTS_WITH(pszVal, "MD_"))
1188 : {
1189 12 : osMetadataFile = CPLFormFilenameSafe(
1190 12 : CPLGetPathSafe(osFilename.c_str()).c_str(),
1191 6 : pszVal + strlen("MD_"), nullptr);
1192 6 : if (VSIStatL(osMetadataFile.c_str(), &sStat) != 0)
1193 : {
1194 6 : osMetadataFile.clear();
1195 : }
1196 : }
1197 : else
1198 : {
1199 0 : osMetadataFile.clear();
1200 : }
1201 : }
1202 : }
1203 : }
1204 6333 : else if (osName != "horizontalCRS" &&
1205 6109 : osName != "horizontalDatumReference" &&
1206 6087 : osName != "horizontalDatumValue" &&
1207 5852 : osName != "productSpecification" &&
1208 5451 : osName != "eastBoundLongitude" &&
1209 5097 : osName != "northBoundLatitude" &&
1210 4743 : osName != "southBoundLatitude" &&
1211 6495 : osName != "westBoundLongitude" && osName != "extentTypeCode" &&
1212 5824 : osName != "verticalCS" && osName != "verticalCoordinateBase" &&
1213 3370 : osName != "verticalDatumReference" &&
1214 4447 : osName != "verticalDatum" && osName != "nameOfHorizontalCRS" &&
1215 3882 : osName != "typeOfHorizontalCRS" && osName != "horizontalCS" &&
1216 2342 : osName != "horizontalDatum" &&
1217 2249 : osName != "nameOfHorizontalDatum" &&
1218 3329 : osName != "primeMeridian" && osName != "spheroid" &&
1219 2136 : osName != "projectionMethod" &&
1220 2008 : osName != "projectionParameter1" &&
1221 1880 : osName != "projectionParameter2" &&
1222 1771 : osName != "projectionParameter3" &&
1223 1701 : osName != "projectionParameter4" &&
1224 1661 : osName != "projectionParameter5" &&
1225 7156 : osName != "falseNorthing" && osName != "falseEasting")
1226 : {
1227 695 : const char *pszVal = poAttr->ReadAsString();
1228 695 : if (pszVal && pszVal[0])
1229 695 : poDS->GDALDataset::SetMetadataItem(osName.c_str(), pszVal);
1230 : }
1231 : }
1232 224 : return osMetadataFile;
1233 : }
1234 :
1235 : /************************************************************************/
1236 : /* S100BaseWriter::S100BaseWriter() */
1237 : /************************************************************************/
1238 :
1239 200 : S100BaseWriter::S100BaseWriter(const char *pszDestFilename,
1240 200 : GDALDataset *poSrcDS, CSLConstList papszOptions)
1241 : : m_osDestFilename(pszDestFilename), m_poSrcDS(poSrcDS),
1242 200 : m_aosOptions(papszOptions)
1243 : {
1244 200 : }
1245 :
1246 : /************************************************************************/
1247 : /* S100BaseWriter::~S100BaseWriter() */
1248 : /************************************************************************/
1249 :
1250 200 : S100BaseWriter::~S100BaseWriter()
1251 : {
1252 : // Check that destructors of derived classes have called themselves their
1253 : // Close implementation
1254 200 : CPLAssert(!m_hdf5);
1255 200 : }
1256 :
1257 : /************************************************************************/
1258 : /* S100BaseWriter::BaseClose() */
1259 : /************************************************************************/
1260 :
1261 284 : bool S100BaseWriter::BaseClose()
1262 : {
1263 284 : bool ret = m_GroupF.clear();
1264 284 : ret = m_valuesGroup.clear() && ret;
1265 284 : ret = m_featureInstanceGroup.clear() && ret;
1266 284 : ret = m_featureGroup.clear() && ret;
1267 284 : ret = m_hdf5.clear() && ret;
1268 284 : return ret;
1269 : }
1270 :
1271 : /************************************************************************/
1272 : /* S100BaseWriter::BaseChecks() */
1273 : /************************************************************************/
1274 :
1275 158 : bool S100BaseWriter::BaseChecks(const char *pszDriverName, bool crsMustBeEPSG,
1276 : bool verticalDatumRequired)
1277 : {
1278 158 : if (m_poSrcDS->GetRasterXSize() < 1 || m_poSrcDS->GetRasterYSize() < 1)
1279 : {
1280 3 : CPLError(CE_Failure, CPLE_NotSupported,
1281 : "Source dataset dimension must be at least 1x1 pixel");
1282 3 : return false;
1283 : }
1284 :
1285 155 : if (m_poSrcDS->GetGeoTransform(m_gt) != CE_None)
1286 : {
1287 6 : CPLError(CE_Failure, CPLE_NotSupported,
1288 : "%s driver requires a source dataset with a geotransform",
1289 : pszDriverName);
1290 6 : return false;
1291 : }
1292 149 : if (m_gt.xrot != 0 || m_gt.yrot != 0)
1293 : {
1294 3 : CPLError(CE_Failure, CPLE_NotSupported,
1295 : "%s driver requires a source dataset with a non-rotated "
1296 : "geotransform",
1297 : pszDriverName);
1298 3 : return false;
1299 : }
1300 :
1301 146 : m_poSRS = m_poSrcDS->GetSpatialRef();
1302 146 : if (!m_poSRS)
1303 : {
1304 3 : CPLError(CE_Failure, CPLE_NotSupported,
1305 : "%s driver requires a source dataset with a CRS",
1306 : pszDriverName);
1307 3 : return false;
1308 : }
1309 :
1310 143 : const char *pszAuthName = m_poSRS->GetAuthorityName();
1311 143 : const char *pszAuthCode = m_poSRS->GetAuthorityCode();
1312 143 : if (pszAuthName && pszAuthCode && EQUAL(pszAuthName, "EPSG"))
1313 : {
1314 79 : m_nEPSGCode = atoi(pszAuthCode);
1315 : }
1316 143 : if (crsMustBeEPSG && m_nEPSGCode == 0)
1317 : {
1318 15 : CPLError(CE_Failure, CPLE_NotSupported,
1319 : "%s driver requires a source dataset whose CRS has an EPSG "
1320 : "identifier",
1321 : pszDriverName);
1322 15 : return false;
1323 : }
1324 :
1325 : const char *pszVerticalDatum =
1326 128 : m_aosOptions.FetchNameValue("VERTICAL_DATUM");
1327 128 : if (!pszVerticalDatum)
1328 : pszVerticalDatum =
1329 43 : m_poSrcDS->GetMetadataItem(S100_VERTICAL_DATUM_EPSG_CODE);
1330 128 : if (!pszVerticalDatum)
1331 42 : pszVerticalDatum = m_poSrcDS->GetMetadataItem(S100_VERTICAL_DATUM_NAME);
1332 128 : if (!pszVerticalDatum)
1333 : {
1334 40 : if (verticalDatumRequired)
1335 : {
1336 4 : CPLError(CE_Failure, CPLE_AppDefined,
1337 : "VERTICAL_DATUM creation option must be specified");
1338 4 : return false;
1339 : }
1340 : }
1341 : else
1342 : {
1343 88 : m_nVerticalDatum =
1344 88 : S100GetVerticalDatumCodeFromNameOrAbbrev(pszVerticalDatum);
1345 88 : if (m_nVerticalDatum <= 0)
1346 : {
1347 4 : auto pjCtxt = OSRGetProjTLSContext();
1348 : PJ *vertical_datum =
1349 4 : proj_create_from_database(pjCtxt, "EPSG", pszVerticalDatum,
1350 : PJ_CATEGORY_DATUM, false, nullptr);
1351 6 : const bool bIsValid = vertical_datum != nullptr &&
1352 2 : proj_get_type(vertical_datum) ==
1353 4 : PJ_TYPE_VERTICAL_REFERENCE_FRAME;
1354 4 : proj_destroy(vertical_datum);
1355 4 : if (bIsValid)
1356 : {
1357 2 : m_nVerticalDatum = atoi(pszVerticalDatum);
1358 : }
1359 : else
1360 : {
1361 2 : CPLError(CE_Failure, CPLE_AppDefined,
1362 : "VERTICAL_DATUM value is invalid");
1363 2 : return false;
1364 : }
1365 : }
1366 : }
1367 :
1368 122 : const std::string osFilename = CPLGetFilename(m_osDestFilename.c_str());
1369 122 : CPLAssert(pszDriverName[0] == 'S');
1370 122 : const char *pszExpectedFilenamePrefix = pszDriverName + 1;
1371 122 : if (!cpl::starts_with(osFilename, pszExpectedFilenamePrefix))
1372 : {
1373 5 : CPLError(CE_Warning, CPLE_AppDefined,
1374 : "%s dataset filenames should start with '%s'", pszDriverName,
1375 : pszExpectedFilenamePrefix);
1376 : }
1377 127 : if (!cpl::ends_with(osFilename, ".h5") &&
1378 5 : !cpl::ends_with(osFilename, ".H5"))
1379 : {
1380 5 : CPLError(CE_Warning, CPLE_AppDefined,
1381 : "%s dataset filenames should have a '.H5' extension",
1382 : pszDriverName);
1383 : }
1384 :
1385 122 : return true;
1386 : }
1387 :
1388 : /************************************************************************/
1389 : /* S100BaseWriter::OpenFileUpdateMode() */
1390 : /************************************************************************/
1391 :
1392 5 : bool S100BaseWriter::OpenFileUpdateMode()
1393 : {
1394 5 : hid_t fapl = H5_CHECK(H5Pcreate(H5P_FILE_ACCESS));
1395 5 : H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr);
1396 5 : m_hdf5.reset(H5Fopen(m_osDestFilename.c_str(), H5F_ACC_RDWR, fapl));
1397 5 : H5Pclose(fapl);
1398 5 : if (!m_hdf5)
1399 : {
1400 0 : CPLError(CE_Failure, CPLE_AppDefined,
1401 : "Cannot open file %s in update mode",
1402 : m_osDestFilename.c_str());
1403 0 : return false;
1404 : }
1405 5 : return true;
1406 : }
1407 :
1408 : /************************************************************************/
1409 : /* S100BaseWriter::CreateFile() */
1410 : /************************************************************************/
1411 :
1412 90 : bool S100BaseWriter::CreateFile()
1413 : {
1414 90 : hid_t fapl = H5_CHECK(H5Pcreate(H5P_FILE_ACCESS));
1415 90 : H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr);
1416 : {
1417 180 : GH5_libhdf5_error_silencer oErrorSilencer;
1418 90 : m_hdf5.reset(H5Fcreate(m_osDestFilename.c_str(), H5F_ACC_TRUNC,
1419 : H5P_DEFAULT, fapl));
1420 : }
1421 90 : H5Pclose(fapl);
1422 90 : if (!m_hdf5)
1423 : {
1424 2 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create file %s",
1425 : m_osDestFilename.c_str());
1426 2 : return false;
1427 : }
1428 88 : return true;
1429 : }
1430 :
1431 : /************************************************************************/
1432 : /* S100BaseWriter::WriteUInt8Value() */
1433 : /************************************************************************/
1434 :
1435 77 : bool S100BaseWriter::WriteUInt8Value(hid_t hGroup, const char *pszName,
1436 : int value)
1437 : {
1438 154 : return GH5_CreateAttribute(hGroup, pszName, H5T_STD_U8LE) &&
1439 154 : GH5_WriteAttribute(hGroup, pszName, value);
1440 : }
1441 :
1442 : /************************************************************************/
1443 : /* S100BaseWriter::WriteUInt16Value() */
1444 : /************************************************************************/
1445 :
1446 4 : bool S100BaseWriter::WriteUInt16Value(hid_t hGroup, const char *pszName,
1447 : int value)
1448 : {
1449 8 : return GH5_CreateAttribute(hGroup, pszName, H5T_STD_U16LE) &&
1450 8 : GH5_WriteAttribute(hGroup, pszName, value);
1451 : }
1452 :
1453 : /************************************************************************/
1454 : /* S100BaseWriter::WriteUInt32Value() */
1455 : /************************************************************************/
1456 :
1457 227 : bool S100BaseWriter::WriteUInt32Value(hid_t hGroup, const char *pszName,
1458 : unsigned value)
1459 : {
1460 454 : return GH5_CreateAttribute(hGroup, pszName, H5T_STD_U32LE) &&
1461 454 : GH5_WriteAttribute(hGroup, pszName, value);
1462 : }
1463 :
1464 : /************************************************************************/
1465 : /* S100BaseWriter::WriteInt32Value() */
1466 : /************************************************************************/
1467 :
1468 203 : bool S100BaseWriter::WriteInt32Value(hid_t hGroup, const char *pszName,
1469 : int value)
1470 : {
1471 406 : return GH5_CreateAttribute(hGroup, pszName, H5T_STD_I32LE) &&
1472 406 : GH5_WriteAttribute(hGroup, pszName, value);
1473 : }
1474 :
1475 : /************************************************************************/
1476 : /* S100BaseWriter::WriteFloat32Value() */
1477 : /************************************************************************/
1478 :
1479 997 : bool S100BaseWriter::WriteFloat32Value(hid_t hGroup, const char *pszName,
1480 : double value)
1481 : {
1482 1994 : return GH5_CreateAttribute(hGroup, pszName, H5T_IEEE_F32LE) &&
1483 1994 : GH5_WriteAttribute(hGroup, pszName, value);
1484 : }
1485 :
1486 : /************************************************************************/
1487 : /* S100BaseWriter::WriteFloat64Value() */
1488 : /************************************************************************/
1489 :
1490 524 : bool S100BaseWriter::WriteFloat64Value(hid_t hGroup, const char *pszName,
1491 : double value)
1492 : {
1493 1048 : return GH5_CreateAttribute(hGroup, pszName, H5T_IEEE_F64LE) &&
1494 1048 : GH5_WriteAttribute(hGroup, pszName, value);
1495 : }
1496 :
1497 : /************************************************************************/
1498 : /* S100BaseWriter::WriteVarLengthStringValue() */
1499 : /************************************************************************/
1500 :
1501 663 : bool S100BaseWriter::WriteVarLengthStringValue(hid_t hGroup,
1502 : const char *pszName,
1503 : const char *pszValue)
1504 : {
1505 1326 : return GH5_CreateAttribute(hGroup, pszName, H5T_C_S1, VARIABLE_LENGTH) &&
1506 1326 : GH5_WriteAttribute(hGroup, pszName, pszValue);
1507 : }
1508 :
1509 : /************************************************************************/
1510 : /* S100BaseWriter::WriteFixedLengthStringValue() */
1511 : /************************************************************************/
1512 :
1513 0 : bool S100BaseWriter::WriteFixedLengthStringValue(hid_t hGroup,
1514 : const char *pszName,
1515 : const char *pszValue)
1516 : {
1517 0 : return GH5_CreateAttribute(hGroup, pszName, H5T_C_S1,
1518 0 : static_cast<unsigned>(strlen(pszValue))) &&
1519 0 : GH5_WriteAttribute(hGroup, pszName, pszValue);
1520 : }
1521 :
1522 : /************************************************************************/
1523 : /* S100BaseWriter::WriteProductSpecification() */
1524 : /************************************************************************/
1525 :
1526 88 : bool S100BaseWriter::WriteProductSpecification(
1527 : const char *pszProductSpecification)
1528 : {
1529 88 : return WriteVarLengthStringValue(m_hdf5, "productSpecification",
1530 88 : pszProductSpecification);
1531 : }
1532 :
1533 : /************************************************************************/
1534 : /* S100BaseWriter::WriteIssueDate() */
1535 : /************************************************************************/
1536 :
1537 88 : bool S100BaseWriter::WriteIssueDate()
1538 : {
1539 88 : const char *pszIssueDate = m_aosOptions.FetchNameValue("ISSUE_DATE");
1540 88 : if (!pszIssueDate)
1541 : {
1542 55 : const char *pszTmp = m_poSrcDS->GetMetadataItem("issueDate");
1543 55 : if (pszTmp && strlen(pszTmp) == 8)
1544 4 : pszIssueDate = pszTmp;
1545 : }
1546 :
1547 176 : std::string osIssueDate; // keep in that scope
1548 88 : if (pszIssueDate)
1549 : {
1550 37 : if (strlen(pszIssueDate) != 8)
1551 0 : CPLError(CE_Warning, CPLE_AppDefined,
1552 : "ISSUE_DATE should be 8 digits: YYYYMMDD");
1553 : }
1554 : else
1555 : {
1556 : time_t now;
1557 51 : time(&now);
1558 : struct tm brokenDown;
1559 51 : CPLUnixTimeToYMDHMS(now, &brokenDown);
1560 51 : osIssueDate = CPLSPrintf("%04d%02d%02d", brokenDown.tm_year + 1900,
1561 51 : brokenDown.tm_mon + 1, brokenDown.tm_mday);
1562 51 : pszIssueDate = osIssueDate.c_str();
1563 : }
1564 :
1565 176 : return WriteVarLengthStringValue(m_hdf5, "issueDate", pszIssueDate);
1566 : }
1567 :
1568 : /************************************************************************/
1569 : /* S100BaseWriter::WriteIssueTime() */
1570 : /************************************************************************/
1571 :
1572 88 : bool S100BaseWriter::WriteIssueTime(bool bAutogenerateFromCurrent)
1573 : {
1574 88 : const char *pszIssueTime = m_aosOptions.FetchNameValue("ISSUE_TIME");
1575 88 : if (!pszIssueTime)
1576 : {
1577 55 : const char *pszTmp = m_poSrcDS->GetMetadataItem("issueTime");
1578 55 : if (pszTmp && strlen(pszTmp) == 7 && pszTmp[6] == 'Z')
1579 4 : pszIssueTime = pszTmp;
1580 : }
1581 88 : std::string osIssueTime; // keep in that scope
1582 88 : if (!pszIssueTime && bAutogenerateFromCurrent)
1583 : {
1584 : time_t now;
1585 36 : time(&now);
1586 : struct tm brokenDown;
1587 36 : CPLUnixTimeToYMDHMS(now, &brokenDown);
1588 : osIssueTime = CPLSPrintf("%02d%02d%02dZ", brokenDown.tm_hour,
1589 36 : brokenDown.tm_min, brokenDown.tm_sec);
1590 36 : pszIssueTime = osIssueTime.c_str();
1591 : }
1592 161 : return !pszIssueTime || pszIssueTime[0] == 0 ||
1593 249 : WriteVarLengthStringValue(m_hdf5, "issueTime", pszIssueTime);
1594 : }
1595 :
1596 : /************************************************************************/
1597 : /* S100BaseWriter::WriteTopLevelBoundingBox() */
1598 : /************************************************************************/
1599 :
1600 85 : bool S100BaseWriter::WriteTopLevelBoundingBox()
1601 : {
1602 :
1603 85 : OGREnvelope sExtent;
1604 85 : if (m_poSrcDS->GetExtentWGS84LongLat(&sExtent) != OGRERR_NONE)
1605 : {
1606 0 : CPLError(CE_Failure, CPLE_AppDefined,
1607 : "Cannot get dataset extent in WGS84 longitude/latitude");
1608 0 : return false;
1609 : }
1610 :
1611 85 : return WriteFloat32Value(m_hdf5, "westBoundLongitude", sExtent.MinX) &&
1612 85 : WriteFloat32Value(m_hdf5, "southBoundLatitude", sExtent.MinY) &&
1613 255 : WriteFloat32Value(m_hdf5, "eastBoundLongitude", sExtent.MaxX) &&
1614 170 : WriteFloat32Value(m_hdf5, "northBoundLatitude", sExtent.MaxY);
1615 : }
1616 :
1617 : /************************************************************************/
1618 : /* S100BaseWriter::WriteHorizontalCRS() */
1619 : /************************************************************************/
1620 :
1621 88 : bool S100BaseWriter::WriteHorizontalCRS()
1622 : {
1623 88 : bool ret = WriteInt32Value(m_hdf5, "horizontalCRS",
1624 88 : m_nEPSGCode > 0 ? m_nEPSGCode : -1);
1625 88 : if (ret && m_nEPSGCode <= 0)
1626 : {
1627 40 : ret = WriteVarLengthStringValue(m_hdf5, "nameOfHorizontalCRS",
1628 40 : m_poSRS->GetName());
1629 : {
1630 80 : GH5_HIDTypeHolder hEnumType(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
1631 40 : ret = ret && hEnumType;
1632 40 : if (ret)
1633 : {
1634 : uint8_t val;
1635 40 : val = 1;
1636 40 : ret = ret && H5_CHECK(H5Tenum_insert(hEnumType, "geodeticCRS2D",
1637 : &val)) >= 0;
1638 40 : val = 2;
1639 40 : ret = ret && H5_CHECK(H5Tenum_insert(hEnumType, "projectedCRS",
1640 : &val)) >= 0;
1641 40 : ret = ret &&
1642 40 : GH5_CreateAttribute(m_hdf5, "typeOfHorizontalCRS",
1643 80 : hEnumType) &&
1644 40 : GH5_WriteAttribute(m_hdf5, "typeOfHorizontalCRS",
1645 40 : m_poSRS->IsGeographic() ? 1 : 2);
1646 : }
1647 : }
1648 :
1649 69 : const int nHorizontalCS = m_poSRS->IsGeographic() ? 6422
1650 29 : : m_poSRS->EPSGTreatsAsNorthingEasting()
1651 29 : ? 4500
1652 40 : : 4400;
1653 40 : ret = ret && WriteInt32Value(m_hdf5, "horizontalCS", nHorizontalCS);
1654 :
1655 : const char *pszDatumKey =
1656 40 : m_poSRS->IsGeographic() ? "GEOGCS|DATUM" : "PROJCS|GEOGCS|DATUM";
1657 40 : const char *pszDatumAuthName = m_poSRS->GetAuthorityName(pszDatumKey);
1658 40 : const char *pszDatumCode = m_poSRS->GetAuthorityCode(pszDatumKey);
1659 34 : const int nDatum = (pszDatumAuthName && pszDatumCode &&
1660 34 : EQUAL(pszDatumAuthName, "EPSG"))
1661 74 : ? atoi(pszDatumCode)
1662 : : -1;
1663 40 : ret = ret && WriteInt32Value(m_hdf5, "horizontalDatum", nDatum);
1664 40 : if (ret && nDatum < 0)
1665 : {
1666 6 : const char *pszDatum = m_poSRS->GetAttrValue(pszDatumKey);
1667 6 : if (!pszDatum)
1668 0 : pszDatum = "unknown";
1669 6 : ret = WriteVarLengthStringValue(m_hdf5, "nameOfHorizontalDatum",
1670 : pszDatum);
1671 :
1672 6 : const char *pszSpheroidKey = m_poSRS->IsGeographic()
1673 6 : ? "GEOGCS|DATUM|SPHEROID"
1674 6 : : "PROJCS|GEOGCS|DATUM|SPHEROID";
1675 : const char *pszSpheroidAuthName =
1676 6 : m_poSRS->GetAuthorityName(pszSpheroidKey);
1677 : const char *pszSpheroidCode =
1678 6 : m_poSRS->GetAuthorityCode(pszSpheroidKey);
1679 6 : const char *pszSpheroidName = m_poSRS->GetAttrValue(pszSpheroidKey);
1680 6 : const int nSpheroid =
1681 3 : (pszSpheroidAuthName && pszSpheroidCode &&
1682 3 : EQUAL(pszSpheroidAuthName, "EPSG"))
1683 12 : ? atoi(pszSpheroidCode)
1684 3 : : (pszSpheroidName && EQUAL(pszSpheroidName, "Bessel 1841"))
1685 6 : ? 7004
1686 : : -1;
1687 6 : if (nSpheroid <= 0)
1688 : {
1689 1 : CPLError(CE_Failure, CPLE_AppDefined,
1690 : "Unknown code for ellipsoid of CRS");
1691 1 : return false;
1692 : }
1693 5 : ret = ret && WriteInt32Value(m_hdf5, "spheroid", nSpheroid);
1694 :
1695 5 : const char *pszPrimeMeridianKey = m_poSRS->IsGeographic()
1696 5 : ? "GEOGCS|PRIMEM"
1697 5 : : "PROJCS|GEOGCS|PRIMEM";
1698 : const char *pszPrimeMeridianAuthName =
1699 5 : m_poSRS->GetAuthorityName(pszPrimeMeridianKey);
1700 : const char *pszPrimeMeridianCode =
1701 5 : m_poSRS->GetAuthorityCode(pszPrimeMeridianKey);
1702 : const char *pszPrimeMeridianName =
1703 5 : m_poSRS->GetAttrValue(pszPrimeMeridianKey);
1704 5 : const int nPrimeMeridian =
1705 2 : (pszPrimeMeridianAuthName && pszPrimeMeridianCode &&
1706 2 : EQUAL(pszPrimeMeridianAuthName, "EPSG"))
1707 10 : ? atoi(pszPrimeMeridianCode)
1708 3 : : (pszPrimeMeridianName && EQUAL(pszPrimeMeridianName, "Ferro"))
1709 6 : ? 8909
1710 : : -1;
1711 5 : if (nPrimeMeridian <= 0)
1712 : {
1713 1 : CPLError(CE_Failure, CPLE_AppDefined,
1714 : "Unknown code for prime meridian of CRS");
1715 1 : return false;
1716 : }
1717 4 : ret =
1718 4 : ret && WriteInt32Value(m_hdf5, "primeMeridian", nPrimeMeridian);
1719 : }
1720 :
1721 38 : const char *pszProjection = m_poSRS->IsProjected()
1722 38 : ? m_poSRS->GetAttrValue("PROJECTION")
1723 38 : : nullptr;
1724 38 : if (pszProjection)
1725 : {
1726 27 : int nProjectionMethod = 0;
1727 27 : double adfParams[] = {std::numeric_limits<double>::quiet_NaN(),
1728 : std::numeric_limits<double>::quiet_NaN(),
1729 : std::numeric_limits<double>::quiet_NaN(),
1730 : std::numeric_limits<double>::quiet_NaN(),
1731 : std::numeric_limits<double>::quiet_NaN()};
1732 27 : if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP))
1733 : {
1734 2 : nProjectionMethod = PROJECTION_METHOD_MERCATOR;
1735 2 : adfParams[0] =
1736 2 : m_poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
1737 2 : adfParams[1] =
1738 2 : m_poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1739 : }
1740 25 : else if (EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
1741 : {
1742 : auto poTmpSRS = std::unique_ptr<OGRSpatialReference>(
1743 2 : m_poSRS->convertToOtherProjection(SRS_PT_MERCATOR_2SP));
1744 2 : nProjectionMethod = PROJECTION_METHOD_MERCATOR;
1745 2 : adfParams[0] =
1746 2 : poTmpSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
1747 2 : adfParams[1] =
1748 2 : poTmpSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1749 : }
1750 23 : else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
1751 : {
1752 2 : nProjectionMethod = PROJECTION_METHOD_TRANSVERSE_MERCATOR;
1753 2 : adfParams[0] =
1754 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
1755 2 : adfParams[1] =
1756 2 : m_poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1757 2 : adfParams[2] =
1758 2 : m_poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
1759 : }
1760 21 : else if (EQUAL(pszProjection,
1761 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
1762 : {
1763 2 : nProjectionMethod = PROJECTION_METHOD_OBLIQUE_MERCATOR;
1764 2 : adfParams[0] =
1765 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0);
1766 2 : adfParams[1] =
1767 2 : m_poSRS->GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
1768 2 : adfParams[2] = m_poSRS->GetNormProjParm(SRS_PP_AZIMUTH, 0.0);
1769 2 : adfParams[3] =
1770 2 : m_poSRS->GetNormProjParm(SRS_PP_RECTIFIED_GRID_ANGLE, 0.0);
1771 2 : adfParams[4] =
1772 2 : m_poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
1773 : }
1774 19 : else if (EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR))
1775 : {
1776 2 : nProjectionMethod = PROJECTION_METHOD_HOTINE_OBLIQUE_MERCATOR;
1777 2 : adfParams[0] =
1778 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0);
1779 2 : adfParams[1] =
1780 2 : m_poSRS->GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
1781 2 : adfParams[2] = m_poSRS->GetNormProjParm(SRS_PP_AZIMUTH, 0.0);
1782 2 : adfParams[3] =
1783 2 : m_poSRS->GetNormProjParm(SRS_PP_RECTIFIED_GRID_ANGLE, 0.0);
1784 2 : adfParams[4] =
1785 2 : m_poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
1786 : }
1787 17 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
1788 : {
1789 2 : nProjectionMethod = PROJECTION_METHOD_LCC_1SP;
1790 2 : adfParams[0] =
1791 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
1792 2 : adfParams[1] =
1793 2 : m_poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1794 2 : adfParams[2] =
1795 2 : m_poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
1796 : }
1797 15 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
1798 : {
1799 2 : nProjectionMethod = PROJECTION_METHOD_LCC_2SP;
1800 2 : adfParams[0] =
1801 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
1802 2 : adfParams[1] =
1803 2 : m_poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1804 2 : adfParams[2] =
1805 2 : m_poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
1806 2 : adfParams[3] =
1807 2 : m_poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0);
1808 : }
1809 13 : else if (EQUAL(pszProjection, SRS_PT_OBLIQUE_STEREOGRAPHIC))
1810 : {
1811 2 : nProjectionMethod = PROJECTION_METHOD_OBLIQUE_STEREOGRAPHIC;
1812 2 : adfParams[0] =
1813 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
1814 2 : adfParams[1] =
1815 2 : m_poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1816 2 : adfParams[2] =
1817 2 : m_poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
1818 : }
1819 11 : else if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
1820 : {
1821 2 : nProjectionMethod = PROJECTION_METHOD_POLAR_STEREOGRAPHIC;
1822 2 : adfParams[0] =
1823 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
1824 2 : adfParams[1] =
1825 2 : m_poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1826 2 : adfParams[2] =
1827 2 : m_poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
1828 : }
1829 9 : else if (EQUAL(pszProjection, SRS_PT_KROVAK))
1830 : {
1831 2 : nProjectionMethod =
1832 : PROJECTION_METHOD_KROVAK_OBLIQUE_CONIC_CONFORMAL;
1833 2 : adfParams[0] =
1834 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
1835 2 : adfParams[1] =
1836 2 : m_poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1837 2 : adfParams[2] = m_poSRS->GetNormProjParm(SRS_PP_AZIMUTH, 0.0);
1838 2 : adfParams[3] =
1839 2 : m_poSRS->GetNormProjParm(SRS_PP_PSEUDO_STD_PARALLEL_1, 0.0);
1840 2 : adfParams[4] =
1841 2 : m_poSRS->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
1842 : }
1843 7 : else if (EQUAL(pszProjection, SRS_PT_POLYCONIC))
1844 : {
1845 2 : nProjectionMethod = PROJECTION_METHOD_AMERICAN_POLYCONIC;
1846 2 : adfParams[0] =
1847 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
1848 2 : adfParams[1] =
1849 2 : m_poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
1850 : }
1851 5 : else if (EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA))
1852 : {
1853 2 : nProjectionMethod = PROJECTION_METHOD_ALBERS_EQUAL_AREA;
1854 2 : adfParams[0] =
1855 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0);
1856 2 : adfParams[1] =
1857 2 : m_poSRS->GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
1858 2 : adfParams[2] =
1859 2 : m_poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
1860 2 : adfParams[3] =
1861 2 : m_poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0);
1862 : }
1863 3 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
1864 : {
1865 2 : nProjectionMethod =
1866 : PROJECTION_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA;
1867 2 : adfParams[0] =
1868 2 : m_poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0);
1869 2 : adfParams[1] =
1870 2 : m_poSRS->GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
1871 : }
1872 : else
1873 : {
1874 1 : CPLError(CE_Failure, CPLE_NotSupported,
1875 : "Projection method %s is not supported by S100",
1876 : pszProjection);
1877 1 : return false;
1878 : }
1879 :
1880 26 : ret = ret && WriteInt32Value(m_hdf5, "projectionMethod",
1881 : nProjectionMethod);
1882 112 : for (int i = 0; i < 5 && !std::isnan(adfParams[i]); ++i)
1883 : {
1884 : const std::string osAttrName =
1885 86 : "projectionParameter" + std::to_string(i + 1);
1886 86 : ret = ret && WriteFloat64Value(m_hdf5, osAttrName.c_str(),
1887 : adfParams[i]);
1888 : }
1889 :
1890 52 : ret = ret && WriteFloat64Value(m_hdf5, "falseNorthing",
1891 26 : m_poSRS->GetNormProjParm(
1892 : SRS_PP_FALSE_NORTHING, 0.0));
1893 52 : ret = ret && WriteFloat64Value(m_hdf5, "falseEasting",
1894 26 : m_poSRS->GetNormProjParm(
1895 : SRS_PP_FALSE_EASTING, 0.0));
1896 : }
1897 : }
1898 85 : return ret;
1899 : }
1900 :
1901 : /************************************************************************/
1902 : /* S100BaseWriter::WriteVerticalCoordinateBase() */
1903 : /************************************************************************/
1904 :
1905 79 : bool S100BaseWriter::WriteVerticalCoordinateBase(int nCode)
1906 : {
1907 79 : GH5_HIDTypeHolder hEnumType(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
1908 79 : bool ret = hEnumType;
1909 79 : if (hEnumType)
1910 : {
1911 : uint8_t val;
1912 79 : val = 1;
1913 79 : ret =
1914 79 : ret && H5_CHECK(H5Tenum_insert(hEnumType, "seaSurface", &val)) >= 0;
1915 79 : val = 2;
1916 158 : ret = ret &&
1917 79 : H5_CHECK(H5Tenum_insert(hEnumType, "verticalDatum", &val)) >= 0;
1918 79 : val = 3;
1919 79 : ret =
1920 79 : ret && H5_CHECK(H5Tenum_insert(hEnumType, "seaBottom", &val)) >= 0;
1921 :
1922 79 : ret =
1923 79 : ret &&
1924 158 : GH5_CreateAttribute(m_hdf5, "verticalCoordinateBase", hEnumType) &&
1925 79 : GH5_WriteAttribute(m_hdf5, "verticalCoordinateBase", nCode);
1926 : }
1927 158 : return ret;
1928 : }
1929 :
1930 : /************************************************************************/
1931 : /* S100BaseWriter::WriteVerticalDatumReference() */
1932 : /************************************************************************/
1933 :
1934 55 : bool S100BaseWriter::WriteVerticalDatumReference(hid_t hGroup, int nCode)
1935 : {
1936 55 : GH5_HIDTypeHolder hEnumType(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
1937 55 : bool ret = hEnumType;
1938 55 : if (hEnumType)
1939 : {
1940 : uint8_t val;
1941 55 : val = 1;
1942 55 : ret = ret && H5_CHECK(H5Tenum_insert(hEnumType, "s100VerticalDatum",
1943 : &val)) >= 0;
1944 55 : val = 2;
1945 55 : ret = ret && H5_CHECK(H5Tenum_insert(hEnumType, "EPSG", &val)) >= 0;
1946 :
1947 55 : ret =
1948 55 : ret &&
1949 110 : GH5_CreateAttribute(hGroup, "verticalDatumReference", hEnumType) &&
1950 55 : GH5_WriteAttribute(hGroup, "verticalDatumReference", nCode);
1951 : }
1952 110 : return ret;
1953 : }
1954 :
1955 : /************************************************************************/
1956 : /* S100BaseWriter::WriteVerticalCS() */
1957 : /************************************************************************/
1958 :
1959 54 : bool S100BaseWriter::WriteVerticalCS(int nCode)
1960 : {
1961 108 : return GH5_CreateAttribute(m_hdf5, "verticalCS", H5T_STD_I32LE) &&
1962 108 : GH5_WriteAttribute(m_hdf5, "verticalCS", nCode);
1963 : }
1964 :
1965 : /************************************************************************/
1966 : /* S100BaseWriter::WriteVerticalDatum() */
1967 : /************************************************************************/
1968 :
1969 56 : bool S100BaseWriter::WriteVerticalDatum(hid_t hGroup, hid_t hType, int nCode)
1970 : {
1971 112 : return GH5_CreateAttribute(hGroup, "verticalDatum", hType) &&
1972 112 : GH5_WriteAttribute(hGroup, "verticalDatum", nCode);
1973 : }
1974 :
1975 : /************************************************************************/
1976 : /* S100BaseWriter::CreateGroupF() */
1977 : /************************************************************************/
1978 :
1979 70 : bool S100BaseWriter::CreateGroupF()
1980 : {
1981 70 : m_GroupF.reset(H5_CHECK(H5Gcreate(m_hdf5, "Group_F", 0)));
1982 70 : return m_GroupF;
1983 : }
1984 :
1985 : /************************************************************************/
1986 : /* S100BaseWriter::CreateFeatureGroup() */
1987 : /************************************************************************/
1988 :
1989 77 : bool S100BaseWriter::CreateFeatureGroup(const char *name)
1990 : {
1991 77 : m_featureGroup.reset(H5_CHECK(H5Gcreate(m_hdf5, name, 0)));
1992 77 : return m_featureGroup;
1993 : }
1994 :
1995 : /************************************************************************/
1996 : /* S100BaseWriter::WriteDataCodingFormat() */
1997 : /************************************************************************/
1998 :
1999 77 : bool S100BaseWriter::WriteDataCodingFormat(hid_t hGroup, int nCode)
2000 : {
2001 77 : GH5_HIDTypeHolder hEnumType(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
2002 77 : bool ret = hEnumType;
2003 77 : if (hEnumType)
2004 : {
2005 77 : uint8_t val = 0;
2006 1386 : for (const char *pszEnumName :
2007 : {"Fixed Stations", "Regular Grid", "Ungeorectified Grid",
2008 : "Moving Platform", "Irregular Grid", "Variable cell size", "TIN",
2009 770 : "Fixed Stations (Stationwise)", "Feature oriented Regular Grid"})
2010 : {
2011 693 : ++val;
2012 1386 : ret = ret &&
2013 693 : H5_CHECK(H5Tenum_insert(hEnumType, pszEnumName, &val)) >= 0;
2014 : }
2015 :
2016 77 : ret = ret &&
2017 154 : GH5_CreateAttribute(hGroup, "dataCodingFormat", hEnumType) &&
2018 77 : GH5_WriteAttribute(hGroup, "dataCodingFormat", nCode);
2019 : }
2020 154 : return ret;
2021 : }
2022 :
2023 : /************************************************************************/
2024 : /* S100BaseWriter::WriteCommonPointRule() */
2025 : /************************************************************************/
2026 :
2027 77 : bool S100BaseWriter::WriteCommonPointRule(hid_t hGroup, int nCode)
2028 : {
2029 77 : GH5_HIDTypeHolder hEnumType(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
2030 77 : bool ret = hEnumType;
2031 77 : if (hEnumType)
2032 : {
2033 77 : uint8_t val = 0;
2034 385 : for (const char *pszEnumName : {"average", "low", "high", "all"})
2035 : {
2036 308 : ++val;
2037 616 : ret = ret &&
2038 308 : H5_CHECK(H5Tenum_insert(hEnumType, pszEnumName, &val)) >= 0;
2039 : }
2040 :
2041 77 : ret = ret &&
2042 154 : GH5_CreateAttribute(hGroup, "commonPointRule", hEnumType) &&
2043 77 : GH5_WriteAttribute(hGroup, "commonPointRule", nCode);
2044 : }
2045 154 : return ret;
2046 : }
2047 :
2048 : /************************************************************************/
2049 : /* S100BaseWriter::WriteDataOffsetCode() */
2050 : /************************************************************************/
2051 :
2052 77 : bool S100BaseWriter::WriteDataOffsetCode(hid_t hGroup, int nCode)
2053 : {
2054 77 : GH5_HIDTypeHolder hEnumType(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
2055 77 : bool ret = hEnumType;
2056 77 : if (hEnumType)
2057 : {
2058 77 : uint8_t val = 0;
2059 770 : for (const char *pszEnumName :
2060 : {"XMin, YMin (\"Lower left\") corner (\"Cell origin\")",
2061 : "XMax, YMax (\"Upper right\") corner",
2062 : "XMax, YMin (\"Lower right\") corner",
2063 : "XMin, YMax (\"Upper left\") corner",
2064 462 : "Barycenter (centroid) of cell"})
2065 : {
2066 385 : ++val;
2067 770 : ret = ret &&
2068 385 : H5_CHECK(H5Tenum_insert(hEnumType, pszEnumName, &val)) >= 0;
2069 : }
2070 :
2071 154 : ret = ret && GH5_CreateAttribute(hGroup, "dataOffsetCode", hEnumType) &&
2072 77 : GH5_WriteAttribute(hGroup, "dataOffsetCode", nCode);
2073 : }
2074 :
2075 154 : return ret;
2076 : }
2077 :
2078 : /************************************************************************/
2079 : /* S100BaseWriter::WriteDimension() */
2080 : /************************************************************************/
2081 :
2082 77 : bool S100BaseWriter::WriteDimension(hid_t hGroup, int nCode)
2083 : {
2084 77 : return WriteUInt8Value(hGroup, "dimension", nCode);
2085 : }
2086 :
2087 : /************************************************************************/
2088 : /* S100BaseWriter::WriteHorizontalPositionUncertainty() */
2089 : /************************************************************************/
2090 :
2091 77 : bool S100BaseWriter::WriteHorizontalPositionUncertainty(hid_t hGroup,
2092 : float fValue)
2093 : {
2094 77 : return WriteFloat32Value(hGroup, "horizontalPositionUncertainty", fValue);
2095 : }
2096 :
2097 : /************************************************************************/
2098 : /* S100BaseWriter::WriteInterpolationType() */
2099 : /************************************************************************/
2100 :
2101 77 : bool S100BaseWriter::WriteInterpolationType(hid_t hGroup, int nCode)
2102 : {
2103 77 : GH5_HIDTypeHolder hEnumType(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
2104 77 : bool ret = hEnumType;
2105 77 : if (hEnumType)
2106 : {
2107 77 : uint8_t val = 0;
2108 77 : constexpr const char *NULL_STRING = nullptr;
2109 770 : for (const char *pszEnumName : {
2110 : "nearestneighbor", // 1
2111 : NULL_STRING, // 2
2112 : NULL_STRING, // 3
2113 : NULL_STRING, // 4
2114 : "bilinear", // 5
2115 : "biquadratic", // 6
2116 : "bicubic", // 7
2117 : NULL_STRING, // 8
2118 : "barycentric", // 9
2119 : "discrete" // 10
2120 847 : })
2121 : {
2122 770 : ++val;
2123 770 : if (pszEnumName)
2124 : {
2125 462 : ret = ret && H5_CHECK(H5Tenum_insert(hEnumType, pszEnumName,
2126 : &val)) >= 0;
2127 : }
2128 : }
2129 :
2130 77 : ret = ret &&
2131 154 : GH5_CreateAttribute(hGroup, "interpolationType", hEnumType) &&
2132 77 : GH5_WriteAttribute(hGroup, "interpolationType", nCode);
2133 : }
2134 154 : return ret;
2135 : }
2136 :
2137 : /************************************************************************/
2138 : /* S100BaseWriter::WriteNumInstances() */
2139 : /************************************************************************/
2140 :
2141 77 : bool S100BaseWriter::WriteNumInstances(hid_t hGroup, hid_t hType,
2142 : int numInstances)
2143 : {
2144 154 : return GH5_CreateAttribute(hGroup, "numInstances", hType) &&
2145 154 : GH5_WriteAttribute(hGroup, "numInstances", numInstances);
2146 : }
2147 :
2148 : /************************************************************************/
2149 : /* S100BaseWriter::WriteSequencingRuleScanDirection() */
2150 : /************************************************************************/
2151 :
2152 77 : bool S100BaseWriter::WriteSequencingRuleScanDirection(hid_t hGroup,
2153 : const char *pszValue)
2154 : {
2155 77 : return WriteVarLengthStringValue(hGroup, "sequencingRule.scanDirection",
2156 77 : pszValue);
2157 : }
2158 :
2159 : /************************************************************************/
2160 : /* S100BaseWriter::WriteSequencingRuleType() */
2161 : /************************************************************************/
2162 :
2163 77 : bool S100BaseWriter::WriteSequencingRuleType(hid_t hGroup, int nCode)
2164 : {
2165 77 : GH5_HIDTypeHolder hEnumType(H5_CHECK(H5Tenum_create(H5T_STD_U8LE)));
2166 77 : bool ret = hEnumType;
2167 77 : if (hEnumType)
2168 : {
2169 77 : uint8_t val = 0;
2170 924 : for (const char *pszEnumName :
2171 : {"linear", "boustrophedonic", "CantorDiagonal", "spiral", "Morton",
2172 539 : "Hilbert"})
2173 : {
2174 462 : ++val;
2175 924 : ret = ret &&
2176 462 : H5_CHECK(H5Tenum_insert(hEnumType, pszEnumName, &val)) >= 0;
2177 : }
2178 :
2179 77 : ret = ret &&
2180 154 : GH5_CreateAttribute(hGroup, "sequencingRule.type", hEnumType) &&
2181 77 : GH5_WriteAttribute(hGroup, "sequencingRule.type", nCode);
2182 : }
2183 154 : return ret;
2184 : }
2185 :
2186 : /************************************************************************/
2187 : /* S100BaseWriter::WriteVerticalUncertainty() */
2188 : /************************************************************************/
2189 :
2190 77 : bool S100BaseWriter::WriteVerticalUncertainty(hid_t hGroup, float fValue)
2191 : {
2192 77 : return WriteFloat32Value(hGroup, "verticalUncertainty", fValue);
2193 : }
2194 :
2195 : /************************************************************************/
2196 : /* S100BaseWriter::WriteOneDimensionalVarLengthStringArray() */
2197 : /************************************************************************/
2198 :
2199 147 : bool S100BaseWriter::WriteOneDimensionalVarLengthStringArray(
2200 : hid_t hGroup, const char *name, CSLConstList values)
2201 : {
2202 147 : bool ret = false;
2203 147 : hsize_t dims[1] = {static_cast<hsize_t>(CSLCount(values))};
2204 294 : GH5_HIDSpaceHolder hSpaceId(H5_CHECK(H5Screate_simple(1, dims, NULL)));
2205 147 : GH5_HIDTypeHolder hTypeId(H5_CHECK(H5Tcopy(H5T_C_S1)));
2206 147 : if (hSpaceId && hTypeId)
2207 : {
2208 294 : ret = H5_CHECK(H5Tset_size(hTypeId, H5T_VARIABLE)) >= 0 &&
2209 147 : H5_CHECK(H5Tset_strpad(hTypeId, H5T_STR_NULLTERM)) >= 0;
2210 294 : GH5_HIDDatasetHolder hDSId;
2211 147 : if (ret)
2212 : {
2213 147 : hDSId.reset(H5_CHECK(
2214 : H5Dcreate(hGroup, name, hTypeId, hSpaceId, H5P_DEFAULT)));
2215 147 : if (hDSId)
2216 147 : ret = H5Dwrite(hDSId, hTypeId, H5S_ALL, H5S_ALL, H5P_DEFAULT,
2217 : values) >= 0;
2218 : }
2219 : }
2220 294 : return ret;
2221 : }
2222 :
2223 : /************************************************************************/
2224 : /* S100BaseWriter::WriteAxisNames() */
2225 : /************************************************************************/
2226 :
2227 77 : bool S100BaseWriter::WriteAxisNames(hid_t hGroup)
2228 : {
2229 77 : const char *axisProjected[] = {"Easting", "Northing", nullptr};
2230 77 : const char *axisGeographic[] = {"Latitude", "Longitude", nullptr};
2231 77 : return WriteOneDimensionalVarLengthStringArray(
2232 : hGroup, "axisNames",
2233 231 : m_poSRS->IsProjected() ? axisProjected : axisGeographic);
2234 : }
2235 :
2236 : /************************************************************************/
2237 : /* S100BaseWriter::CreateFeatureInstanceGroup() */
2238 : /************************************************************************/
2239 :
2240 82 : bool S100BaseWriter::CreateFeatureInstanceGroup(const char *name)
2241 : {
2242 82 : CPLAssert(m_featureGroup);
2243 82 : m_featureInstanceGroup.reset(H5_CHECK(H5Gcreate(m_featureGroup, name, 0)));
2244 82 : return m_featureInstanceGroup;
2245 : }
2246 :
2247 : /************************************************************************/
2248 : /* S100BaseWriter::WriteFIGGridRelatedParameters() */
2249 : /************************************************************************/
2250 :
2251 82 : bool S100BaseWriter::WriteFIGGridRelatedParameters(hid_t hGroup)
2252 : {
2253 : // From pixel-corner convention to pixel-center convention
2254 82 : const double dfMinX = m_gt.xorig + m_gt.xscale / 2;
2255 82 : const double dfMinY = m_gt.yscale < 0
2256 82 : ? m_gt.yorig +
2257 14 : m_gt.yscale * m_poSrcDS->GetRasterYSize() -
2258 14 : m_gt.yscale / 2
2259 68 : : m_gt.yorig + m_gt.yscale / 2;
2260 : const double dfMaxX =
2261 82 : dfMinX + (m_poSrcDS->GetRasterXSize() - 1) * m_gt.xscale;
2262 : const double dfMaxY =
2263 82 : dfMinY + (m_poSrcDS->GetRasterYSize() - 1) * std::fabs(m_gt.yscale);
2264 :
2265 82 : return WriteFloat32Value(hGroup, "westBoundLongitude", dfMinX) &&
2266 82 : WriteFloat32Value(hGroup, "southBoundLatitude", dfMinY) &&
2267 82 : WriteFloat32Value(hGroup, "eastBoundLongitude", dfMaxX) &&
2268 82 : WriteFloat32Value(hGroup, "northBoundLatitude", dfMaxY) &&
2269 82 : WriteFloat64Value(hGroup, "gridOriginLongitude", dfMinX) &&
2270 82 : WriteFloat64Value(hGroup, "gridOriginLatitude", dfMinY) &&
2271 82 : WriteFloat64Value(hGroup, "gridSpacingLongitudinal", m_gt.xscale) &&
2272 82 : WriteFloat64Value(hGroup, "gridSpacingLatitudinal",
2273 82 : std::fabs(m_gt.yscale)) &&
2274 82 : WriteUInt32Value(hGroup, "numPointsLongitudinal",
2275 82 : m_poSrcDS->GetRasterXSize()) &&
2276 82 : WriteUInt32Value(hGroup, "numPointsLatitudinal",
2277 246 : m_poSrcDS->GetRasterYSize()) &&
2278 164 : WriteVarLengthStringValue(hGroup, "startSequence", "0,0");
2279 : }
2280 :
2281 : /************************************************************************/
2282 : /* S100BaseWriter::WriteNumGRP() */
2283 : /************************************************************************/
2284 :
2285 82 : bool S100BaseWriter::WriteNumGRP(hid_t hGroup, hid_t hType, int numGRP)
2286 : {
2287 164 : return GH5_CreateAttribute(hGroup, "numGRP", hType) &&
2288 164 : GH5_WriteAttribute(hGroup, "numGRP", numGRP);
2289 : }
2290 :
2291 : /************************************************************************/
2292 : /* S100BaseWriter::CreateValuesGroup() */
2293 : /************************************************************************/
2294 :
2295 86 : bool S100BaseWriter::CreateValuesGroup(const char *name)
2296 : {
2297 86 : CPLAssert(m_featureInstanceGroup);
2298 86 : m_valuesGroup.reset(H5_CHECK(H5Gcreate(m_featureInstanceGroup, name, 0)));
2299 86 : return m_valuesGroup;
2300 : }
2301 :
2302 : /************************************************************************/
2303 : /* S100BaseWriter::WriteGroupFDataset() */
2304 : /************************************************************************/
2305 :
2306 85 : bool S100BaseWriter::WriteGroupFDataset(
2307 : const char *name,
2308 : const std::vector<std::array<const char *, GROUP_F_DATASET_FIELD_COUNT>>
2309 : &rows)
2310 : {
2311 : GH5_HIDTypeHolder hDataType(H5_CHECK(
2312 170 : H5Tcreate(H5T_COMPOUND, GROUP_F_DATASET_FIELD_COUNT * sizeof(char *))));
2313 170 : GH5_HIDTypeHolder hVarLengthType(H5_CHECK(H5Tcopy(H5T_C_S1)));
2314 : bool bRet =
2315 170 : hDataType && hVarLengthType &&
2316 85 : H5_CHECK(H5Tset_size(hVarLengthType, H5T_VARIABLE)) >= 0 &&
2317 85 : H5_CHECK(H5Tset_strpad(hVarLengthType, H5T_STR_NULLTERM)) >= 0 &&
2318 85 : H5_CHECK(H5Tinsert(hDataType, "code", 0 * sizeof(char *),
2319 85 : hVarLengthType)) >= 0 &&
2320 85 : H5_CHECK(H5Tinsert(hDataType, "name", 1 * sizeof(char *),
2321 85 : hVarLengthType)) >= 0 &&
2322 85 : H5_CHECK(H5Tinsert(hDataType, "uom.name", 2 * sizeof(char *),
2323 85 : hVarLengthType)) >= 0 &&
2324 85 : H5_CHECK(H5Tinsert(hDataType, "fillValue", 3 * sizeof(char *),
2325 85 : hVarLengthType)) >= 0 &&
2326 85 : H5_CHECK(H5Tinsert(hDataType, "datatype", 4 * sizeof(char *),
2327 85 : hVarLengthType)) >= 0 &&
2328 85 : H5_CHECK(H5Tinsert(hDataType, "lower", 5 * sizeof(char *),
2329 85 : hVarLengthType)) >= 0 &&
2330 85 : H5_CHECK(H5Tinsert(hDataType, "upper", 6 * sizeof(char *),
2331 170 : hVarLengthType)) >= 0 &&
2332 85 : H5_CHECK(H5Tinsert(hDataType, "closure", 7 * sizeof(char *),
2333 85 : hVarLengthType)) >= 0;
2334 :
2335 85 : hsize_t dims[] = {static_cast<hsize_t>(rows.size())};
2336 170 : GH5_HIDSpaceHolder hDataSpace(H5_CHECK(H5Screate_simple(1, dims, nullptr)));
2337 85 : bRet = bRet && hDataSpace;
2338 170 : GH5_HIDDatasetHolder hDatasetID;
2339 85 : if (bRet)
2340 : {
2341 85 : hDatasetID.reset(H5_CHECK(
2342 : H5Dcreate(m_GroupF, name, hDataType, hDataSpace, H5P_DEFAULT)));
2343 85 : bRet = hDatasetID;
2344 : }
2345 170 : GH5_HIDSpaceHolder hFileSpace;
2346 85 : if (bRet)
2347 : {
2348 85 : hFileSpace.reset(H5_CHECK(H5Dget_space(hDatasetID)));
2349 85 : bRet = hFileSpace;
2350 : }
2351 :
2352 85 : hsize_t count[] = {1};
2353 85 : GH5_HIDSpaceHolder hMemSpace(H5_CHECK(H5Screate_simple(1, count, nullptr)));
2354 85 : bRet = bRet && hMemSpace;
2355 :
2356 85 : H5OFFSET_TYPE nOffset = 0;
2357 239 : for (const auto &row : rows)
2358 : {
2359 154 : H5OFFSET_TYPE offset[] = {nOffset};
2360 308 : bRet = bRet &&
2361 154 : H5_CHECK(H5Sselect_hyperslab(hFileSpace, H5S_SELECT_SET, offset,
2362 308 : nullptr, count, nullptr)) >= 0 &&
2363 154 : H5_CHECK(H5Dwrite(hDatasetID, hDataType, hMemSpace, hFileSpace,
2364 : H5P_DEFAULT, row.data())) >= 0;
2365 154 : ++nOffset;
2366 : }
2367 :
2368 170 : return bRet;
2369 : }
|
