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