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 "s100.h"
14 : #include "hdf5dataset.h"
15 :
16 : #include "proj.h"
17 : #include "proj_experimental.h"
18 : #include "proj_constants.h"
19 : #include "ogr_proj_p.h"
20 :
21 : #include <algorithm>
22 : #include <cmath>
23 :
24 : /************************************************************************/
25 : /* S100BaseDataset() */
26 : /************************************************************************/
27 :
28 55 : S100BaseDataset::S100BaseDataset(const std::string &osFilename)
29 55 : : m_osFilename(osFilename)
30 : {
31 55 : }
32 :
33 : /************************************************************************/
34 : /* Init() */
35 : /************************************************************************/
36 :
37 55 : bool S100BaseDataset::Init()
38 : {
39 : // Open the file as an HDF5 file.
40 55 : hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
41 55 : H5Pset_driver(fapl, HDF5GetFileDriver(), nullptr);
42 55 : hid_t hHDF5 = H5Fopen(m_osFilename.c_str(), H5F_ACC_RDONLY, fapl);
43 55 : H5Pclose(fapl);
44 55 : if (hHDF5 < 0)
45 0 : return false;
46 :
47 110 : auto poSharedResources = GDAL::HDF5SharedResources::Create(m_osFilename);
48 55 : poSharedResources->m_hHDF5 = hHDF5;
49 :
50 55 : m_poRootGroup = HDF5Dataset::OpenGroup(poSharedResources);
51 55 : if (m_poRootGroup == nullptr)
52 0 : return false;
53 :
54 55 : S100ReadSRS(m_poRootGroup.get(), m_oSRS);
55 :
56 55 : S100ReadVerticalDatum(this, m_poRootGroup.get());
57 :
58 : m_osMetadataFile =
59 55 : S100ReadMetadata(this, m_osFilename, m_poRootGroup.get());
60 :
61 55 : return true;
62 : }
63 :
64 : /************************************************************************/
65 : /* GetGeoTransform() */
66 : /************************************************************************/
67 :
68 12 : CPLErr S100BaseDataset::GetGeoTransform(GDALGeoTransform >) const
69 :
70 : {
71 12 : if (m_bHasGT)
72 : {
73 12 : gt = m_gt;
74 12 : return CE_None;
75 : }
76 :
77 0 : return GDALPamDataset::GetGeoTransform(gt);
78 : }
79 :
80 : /************************************************************************/
81 : /* GetSpatialRef() */
82 : /************************************************************************/
83 :
84 24 : const OGRSpatialReference *S100BaseDataset::GetSpatialRef() const
85 : {
86 24 : if (!m_oSRS.IsEmpty())
87 24 : return &m_oSRS;
88 0 : return GDALPamDataset::GetSpatialRef();
89 : }
90 :
91 : /************************************************************************/
92 : /* GetFileList() */
93 : /************************************************************************/
94 :
95 4 : char **S100BaseDataset::GetFileList()
96 : {
97 4 : char **papszFileList = GDALPamDataset::GetFileList();
98 4 : if (!m_osMetadataFile.empty())
99 4 : papszFileList = CSLAddString(papszFileList, m_osMetadataFile.c_str());
100 4 : return papszFileList;
101 : }
102 :
103 : /************************************************************************/
104 : /* S100ReadSRS() */
105 : /************************************************************************/
106 :
107 85 : bool S100ReadSRS(const GDALGroup *poRootGroup, OGRSpatialReference &oSRS)
108 : {
109 85 : constexpr int PROJECTION_METHOD_MERCATOR = 9805;
110 : static_assert(PROJECTION_METHOD_MERCATOR ==
111 : EPSG_CODE_METHOD_MERCATOR_VARIANT_B);
112 85 : constexpr int PROJECTION_METHOD_TRANSVERSE_MERCATOR = 9807;
113 : static_assert(PROJECTION_METHOD_TRANSVERSE_MERCATOR ==
114 : EPSG_CODE_METHOD_TRANSVERSE_MERCATOR);
115 85 : constexpr int PROJECTION_METHOD_OBLIQUE_MERCATOR = 9815;
116 : static_assert(PROJECTION_METHOD_OBLIQUE_MERCATOR ==
117 : EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_B);
118 85 : constexpr int PROJECTION_METHOD_HOTINE_OBLIQUE_MERCATOR = 9812;
119 : static_assert(PROJECTION_METHOD_HOTINE_OBLIQUE_MERCATOR ==
120 : EPSG_CODE_METHOD_HOTINE_OBLIQUE_MERCATOR_VARIANT_A);
121 85 : constexpr int PROJECTION_METHOD_LCC_1SP = 9801;
122 : static_assert(PROJECTION_METHOD_LCC_1SP ==
123 : EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_1SP);
124 85 : constexpr int PROJECTION_METHOD_LCC_2SP = 9802;
125 : static_assert(PROJECTION_METHOD_LCC_2SP ==
126 : EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_2SP);
127 85 : constexpr int PROJECTION_METHOD_OBLIQUE_STEREOGRAPHIC = 9809;
128 : static_assert(PROJECTION_METHOD_OBLIQUE_STEREOGRAPHIC ==
129 : EPSG_CODE_METHOD_OBLIQUE_STEREOGRAPHIC);
130 85 : constexpr int PROJECTION_METHOD_POLAR_STEREOGRAPHIC = 9810;
131 : static_assert(PROJECTION_METHOD_POLAR_STEREOGRAPHIC ==
132 : EPSG_CODE_METHOD_POLAR_STEREOGRAPHIC_VARIANT_A);
133 85 : constexpr int PROJECTION_METHOD_KROVAK_OBLIQUE_CONIC_CONFORMAL = 9819;
134 : static_assert(PROJECTION_METHOD_KROVAK_OBLIQUE_CONIC_CONFORMAL ==
135 : EPSG_CODE_METHOD_KROVAK);
136 85 : constexpr int PROJECTION_METHOD_AMERICAN_POLYCONIC = 9818;
137 : static_assert(PROJECTION_METHOD_AMERICAN_POLYCONIC ==
138 : EPSG_CODE_METHOD_AMERICAN_POLYCONIC);
139 85 : constexpr int PROJECTION_METHOD_ALBERS_EQUAL_AREA = 9822;
140 : static_assert(PROJECTION_METHOD_ALBERS_EQUAL_AREA ==
141 : EPSG_CODE_METHOD_ALBERS_EQUAL_AREA);
142 85 : constexpr int PROJECTION_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA = 9820;
143 : static_assert(PROJECTION_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA ==
144 : EPSG_CODE_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA);
145 :
146 : // Get SRS
147 85 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
148 255 : auto poHorizontalCRS = poRootGroup->GetAttribute("horizontalCRS");
149 147 : if (poHorizontalCRS &&
150 147 : poHorizontalCRS->GetDataType().GetClass() == GEDTC_NUMERIC)
151 : {
152 : // horizontalCRS is v2.2
153 62 : const int nHorizontalCRS = poHorizontalCRS->ReadAsInt();
154 62 : if (nHorizontalCRS > 0)
155 : {
156 48 : if (oSRS.importFromEPSG(nHorizontalCRS) != OGRERR_NONE)
157 : {
158 0 : oSRS.Clear();
159 : }
160 : }
161 : else
162 : {
163 : auto poNameOfHorizontalCRS =
164 28 : poRootGroup->GetAttribute("nameOfHorizontalCRS");
165 : auto poTypeOfHorizontalCRS =
166 28 : poRootGroup->GetAttribute("typeOfHorizontalCRS");
167 28 : auto poHorizontalCS = poRootGroup->GetAttribute("horizontalCS");
168 : auto poHorizontalDatum =
169 28 : poRootGroup->GetAttribute("horizontalDatum");
170 14 : if (poNameOfHorizontalCRS &&
171 14 : poNameOfHorizontalCRS->GetDataType().GetClass() ==
172 28 : GEDTC_STRING &&
173 28 : poTypeOfHorizontalCRS &&
174 14 : poTypeOfHorizontalCRS->GetDataType().GetClass() ==
175 14 : GEDTC_NUMERIC &&
176 28 : poHorizontalCS &&
177 28 : poHorizontalCS->GetDataType().GetClass() == GEDTC_NUMERIC &&
178 42 : poHorizontalDatum &&
179 14 : poHorizontalDatum->GetDataType().GetClass() == GEDTC_NUMERIC)
180 : {
181 14 : const int nCRSType = poTypeOfHorizontalCRS->ReadAsInt();
182 14 : constexpr int CRS_TYPE_GEOGRAPHIC = 1;
183 14 : constexpr int CRS_TYPE_PROJECTED = 2;
184 :
185 14 : PJ *crs = nullptr;
186 14 : const char *pszCRSName = poNameOfHorizontalCRS->ReadAsString();
187 :
188 14 : const int nHorizontalCS = poHorizontalCS->ReadAsInt();
189 14 : const int nHorizontalDatum = poHorizontalDatum->ReadAsInt();
190 :
191 14 : auto pjCtxt = OSRGetProjTLSContext();
192 14 : if (nHorizontalDatum < 0)
193 : {
194 : auto poNameOfHorizontalDatum =
195 2 : poRootGroup->GetAttribute("nameOfHorizontalDatum");
196 : auto poPrimeMeridian =
197 2 : poRootGroup->GetAttribute("primeMeridian");
198 2 : auto poSpheroid = poRootGroup->GetAttribute("spheroid");
199 1 : if (poNameOfHorizontalDatum &&
200 1 : poNameOfHorizontalDatum->GetDataType().GetClass() ==
201 2 : GEDTC_STRING &&
202 2 : poPrimeMeridian &&
203 1 : poPrimeMeridian->GetDataType().GetClass() ==
204 1 : GEDTC_NUMERIC &&
205 3 : poSpheroid &&
206 1 : poSpheroid->GetDataType().GetClass() == GEDTC_NUMERIC)
207 : {
208 : const char *pszDatumName =
209 1 : poNameOfHorizontalDatum->ReadAsString();
210 1 : const char *pszGeogCRSName =
211 : nCRSType == CRS_TYPE_PROJECTED
212 2 : ? (pszDatumName ? pszDatumName : "unknown")
213 1 : : (pszCRSName ? pszCRSName : "unknown");
214 :
215 1 : const int nPrimeMeridian = poPrimeMeridian->ReadAsInt();
216 1 : const int nSpheroid = poSpheroid->ReadAsInt();
217 1 : PJ *prime_meridian = proj_create_from_database(
218 : pjCtxt, "EPSG", CPLSPrintf("%d", nPrimeMeridian),
219 : PJ_CATEGORY_PRIME_MERIDIAN, false, nullptr);
220 1 : if (!prime_meridian)
221 0 : return false;
222 1 : PJ *spheroid = proj_create_from_database(
223 : pjCtxt, "EPSG", CPLSPrintf("%d", nSpheroid),
224 : PJ_CATEGORY_ELLIPSOID, false, nullptr);
225 1 : if (!spheroid)
226 : {
227 0 : proj_destroy(prime_meridian);
228 0 : return false;
229 : }
230 :
231 1 : double semiMajorMetre = 0;
232 1 : double invFlattening = 0;
233 1 : proj_ellipsoid_get_parameters(pjCtxt, spheroid,
234 : &semiMajorMetre, nullptr,
235 : nullptr, &invFlattening);
236 :
237 1 : double primeMeridianOffset = 0;
238 1 : double primeMeridianUnitConv = 1;
239 1 : const char *primeMeridianUnitName = nullptr;
240 1 : proj_prime_meridian_get_parameters(
241 : pjCtxt, prime_meridian, &primeMeridianOffset,
242 : &primeMeridianUnitConv, &primeMeridianUnitName);
243 :
244 1 : PJ *ellipsoidalCS = proj_create_ellipsoidal_2D_cs(
245 : pjCtxt, PJ_ELLPS2D_LATITUDE_LONGITUDE, nullptr,
246 : 0.0);
247 :
248 1 : crs = proj_create_geographic_crs(
249 : pjCtxt, pszGeogCRSName,
250 : pszDatumName ? pszDatumName : "unknown",
251 : proj_get_name(spheroid), semiMajorMetre,
252 : invFlattening, proj_get_name(prime_meridian),
253 : primeMeridianOffset, primeMeridianUnitName,
254 : primeMeridianUnitConv, ellipsoidalCS);
255 :
256 1 : proj_destroy(ellipsoidalCS);
257 1 : proj_destroy(spheroid);
258 1 : proj_destroy(prime_meridian);
259 : }
260 : else
261 : {
262 0 : CPLError(CE_Warning, CPLE_AppDefined,
263 : "Cannot establish CRS because "
264 : "nameOfHorizontalDatum, primeMeridian and/or "
265 : "spheroid are missing");
266 0 : return false;
267 : }
268 : }
269 : else
270 : {
271 13 : PJ *datum = proj_create_from_database(
272 : pjCtxt, "EPSG", CPLSPrintf("%d", nHorizontalDatum),
273 : PJ_CATEGORY_DATUM, false, nullptr);
274 13 : if (!datum)
275 0 : return false;
276 13 : const char *pszDatumName = proj_get_name(datum);
277 13 : const char *pszGeogCRSName =
278 : nCRSType == CRS_TYPE_PROJECTED
279 14 : ? pszDatumName
280 1 : : (pszCRSName ? pszCRSName : "unknown");
281 :
282 13 : PJ *ellipsoidalCS = proj_create_ellipsoidal_2D_cs(
283 : pjCtxt, PJ_ELLPS2D_LATITUDE_LONGITUDE, nullptr, 0.0);
284 :
285 13 : crs = proj_create_geographic_crs_from_datum(
286 : pjCtxt, pszGeogCRSName, datum, ellipsoidalCS);
287 :
288 13 : proj_destroy(ellipsoidalCS);
289 13 : proj_destroy(datum);
290 : }
291 :
292 14 : if (nCRSType == CRS_TYPE_PROJECTED)
293 : {
294 : auto poProjectionMethod =
295 24 : poRootGroup->GetAttribute("projectionMethod");
296 24 : if (poProjectionMethod &&
297 12 : poProjectionMethod->GetDataType().GetClass() ==
298 12 : GEDTC_NUMERIC)
299 : {
300 : const int nProjectionMethod =
301 12 : poProjectionMethod->ReadAsInt();
302 :
303 : auto poFalseEasting =
304 24 : poRootGroup->GetAttribute("falseEasting");
305 : const double falseEasting =
306 12 : poFalseEasting &&
307 12 : poFalseEasting->GetDataType().GetClass() ==
308 12 : GEDTC_NUMERIC
309 24 : ? poFalseEasting->ReadAsDouble()
310 12 : : 0.0;
311 :
312 : auto poFalseNorthing =
313 24 : poRootGroup->GetAttribute("falseNorthing");
314 : const double falseNorthing =
315 12 : poFalseNorthing &&
316 12 : poFalseNorthing->GetDataType().GetClass() ==
317 12 : GEDTC_NUMERIC
318 24 : ? poFalseNorthing->ReadAsDouble()
319 12 : : 0.0;
320 :
321 41 : const auto getParameter = [poRootGroup](int nParam)
322 : {
323 : auto poParam = poRootGroup->GetAttribute(
324 123 : "projectionParameter" + std::to_string(nParam));
325 41 : if (poParam && poParam->GetDataType().GetClass() ==
326 41 : GEDTC_NUMERIC)
327 41 : return poParam->ReadAsDouble();
328 0 : CPLError(
329 : CE_Warning, CPLE_AppDefined, "%s",
330 0 : std::string(
331 : "Missing attribute projectionParameter")
332 0 : .append(std::to_string(nParam))
333 : .c_str());
334 0 : return 0.0;
335 12 : };
336 :
337 12 : PJ *conv = nullptr;
338 12 : switch (nProjectionMethod)
339 : {
340 1 : case PROJECTION_METHOD_MERCATOR:
341 : {
342 : conv =
343 1 : proj_create_conversion_mercator_variant_b(
344 : pjCtxt,
345 : /* latitude_first_parallel = */
346 : getParameter(1),
347 : /* center_long = */ getParameter(2),
348 : falseEasting, falseNorthing, nullptr,
349 : 0.0, nullptr, 0.0);
350 1 : break;
351 : }
352 :
353 1 : case PROJECTION_METHOD_TRANSVERSE_MERCATOR:
354 : {
355 : conv =
356 1 : proj_create_conversion_transverse_mercator(
357 : pjCtxt,
358 : /* center_lat = */ getParameter(1),
359 : /* center_long = */ getParameter(2),
360 : /* scale = */ getParameter(3),
361 : falseEasting, falseNorthing, nullptr,
362 : 0.0, nullptr, 0.0);
363 1 : break;
364 : }
365 :
366 1 : case PROJECTION_METHOD_OBLIQUE_MERCATOR:
367 : {
368 : conv =
369 1 : proj_create_conversion_hotine_oblique_mercator_variant_b(
370 : pjCtxt,
371 : /* latitude_projection_centre = */
372 : getParameter(1),
373 : /* longitude_projection_centre = */
374 : getParameter(2),
375 : /* azimuth_initial_line = */
376 : getParameter(3),
377 : /* angle_from_rectified_to_skrew_grid = */
378 : getParameter(4),
379 : /* scale = */ getParameter(5),
380 : falseEasting, falseNorthing, nullptr,
381 : 0.0, nullptr, 0.0);
382 1 : break;
383 : }
384 :
385 1 : case PROJECTION_METHOD_HOTINE_OBLIQUE_MERCATOR:
386 : {
387 : conv =
388 1 : proj_create_conversion_hotine_oblique_mercator_variant_a(
389 : pjCtxt,
390 : /* latitude_projection_centre = */
391 : getParameter(1),
392 : /* longitude_projection_centre = */
393 : getParameter(2),
394 : /* azimuth_initial_line = */
395 : getParameter(3),
396 : /* angle_from_rectified_to_skrew_grid = */
397 : getParameter(4),
398 : /* scale = */ getParameter(5),
399 : falseEasting, falseNorthing, nullptr,
400 : 0.0, nullptr, 0.0);
401 1 : break;
402 : }
403 :
404 1 : case PROJECTION_METHOD_LCC_1SP:
405 : {
406 : conv =
407 1 : proj_create_conversion_lambert_conic_conformal_1sp(
408 : pjCtxt,
409 : /* center_lat = */ getParameter(1),
410 : /* center_long = */ getParameter(2),
411 : /* scale = */ getParameter(3),
412 : falseEasting, falseNorthing, nullptr,
413 : 0.0, nullptr, 0.0);
414 1 : break;
415 : }
416 :
417 1 : case PROJECTION_METHOD_LCC_2SP:
418 : {
419 : conv =
420 1 : proj_create_conversion_lambert_conic_conformal_2sp(
421 : pjCtxt,
422 : /* latitude_false_origin = */
423 : getParameter(1),
424 : /* longitude_false_origin = */
425 : getParameter(2),
426 : /* latitude_first_parallel = */
427 : getParameter(3),
428 : /* latitude_second_parallel = */
429 : getParameter(4), falseEasting,
430 : falseNorthing, nullptr, 0.0, nullptr,
431 : 0.0);
432 1 : break;
433 : }
434 :
435 1 : case PROJECTION_METHOD_OBLIQUE_STEREOGRAPHIC:
436 : {
437 : conv =
438 1 : proj_create_conversion_oblique_stereographic(
439 : pjCtxt,
440 : /* center_lat = */ getParameter(1),
441 : /* center_long = */ getParameter(2),
442 : /* scale = */ getParameter(3),
443 : falseEasting, falseNorthing, nullptr,
444 : 0.0, nullptr, 0.0);
445 1 : break;
446 : }
447 :
448 1 : case PROJECTION_METHOD_POLAR_STEREOGRAPHIC:
449 : {
450 : conv =
451 1 : proj_create_conversion_polar_stereographic_variant_a(
452 : pjCtxt,
453 : /* center_lat = */ getParameter(1),
454 : /* center_long = */ getParameter(2),
455 : /* scale = */ getParameter(3),
456 : falseEasting, falseNorthing, nullptr,
457 : 0.0, nullptr, 0.0);
458 1 : break;
459 : }
460 :
461 1 : case PROJECTION_METHOD_KROVAK_OBLIQUE_CONIC_CONFORMAL:
462 : {
463 1 : conv = proj_create_conversion_krovak(
464 : pjCtxt,
465 : /* latitude_projection_centre = */
466 : getParameter(1),
467 : /* longitude_of_origin = */ getParameter(2),
468 : /* colatitude_cone_axis = */
469 : getParameter(3),
470 : /* latitude_pseudo_standard_parallel = */
471 : getParameter(4),
472 : /* scale_factor_pseudo_standard_parallel = */
473 : getParameter(5), falseEasting,
474 : falseNorthing, nullptr, 0.0, nullptr, 0.0);
475 1 : break;
476 : }
477 :
478 1 : case PROJECTION_METHOD_AMERICAN_POLYCONIC:
479 : {
480 : conv =
481 1 : proj_create_conversion_american_polyconic(
482 : pjCtxt,
483 : /* center_lat = */ getParameter(1),
484 : /* center_long = */ getParameter(2),
485 : falseEasting, falseNorthing, nullptr,
486 : 0.0, nullptr, 0.0);
487 1 : break;
488 : }
489 :
490 1 : case PROJECTION_METHOD_ALBERS_EQUAL_AREA:
491 : {
492 1 : conv = proj_create_conversion_albers_equal_area(
493 : pjCtxt,
494 : /* latitude_false_origin = */
495 : getParameter(1),
496 : /* longitude_false_origin = */
497 : getParameter(2),
498 : /* latitude_first_parallel = */
499 : getParameter(3),
500 : /* latitude_second_parallel = */
501 : getParameter(4), falseEasting,
502 : falseNorthing, nullptr, 0.0, nullptr, 0.0);
503 1 : break;
504 : }
505 :
506 1 : case PROJECTION_METHOD_LAMBERT_AZIMUTHAL_EQUAL_AREA:
507 : {
508 : conv =
509 1 : proj_create_conversion_lambert_azimuthal_equal_area(
510 : pjCtxt,
511 : /* latitude_nat_origin = */
512 : getParameter(1),
513 : /* longitude_nat_origin = */
514 : getParameter(2), falseEasting,
515 : falseNorthing, nullptr, 0.0, nullptr,
516 : 0.0);
517 1 : break;
518 : }
519 :
520 0 : default:
521 : {
522 0 : CPLError(CE_Warning, CPLE_AppDefined,
523 : "Cannot establish CRS because "
524 : "projectionMethod = %d is unhandled",
525 : nProjectionMethod);
526 0 : proj_destroy(crs);
527 0 : return false;
528 : }
529 : }
530 :
531 12 : constexpr int HORIZONTAL_CS_EASTING_NORTHING_METRE =
532 : 4400;
533 : // constexpr int HORIZONTAL_CS_NORTHING_EASTING_METRE = 4500;
534 24 : PJ *cartesianCS = proj_create_cartesian_2D_cs(
535 : pjCtxt,
536 : nHorizontalCS ==
537 : HORIZONTAL_CS_EASTING_NORTHING_METRE
538 12 : ? PJ_CART2D_EASTING_NORTHING
539 : : PJ_CART2D_NORTHING_EASTING,
540 : nullptr, 0.0);
541 :
542 12 : PJ *projCRS = proj_create_projected_crs(
543 : pjCtxt, pszCRSName ? pszCRSName : "unknown", crs,
544 : conv, cartesianCS);
545 12 : proj_destroy(crs);
546 12 : crs = projCRS;
547 12 : proj_destroy(conv);
548 12 : proj_destroy(cartesianCS);
549 : }
550 : else
551 : {
552 0 : CPLError(CE_Warning, CPLE_AppDefined,
553 : "Cannot establish CRS because "
554 : "projectionMethod is missing");
555 0 : proj_destroy(crs);
556 0 : return false;
557 : }
558 : }
559 2 : else if (nCRSType != CRS_TYPE_GEOGRAPHIC)
560 : {
561 0 : CPLError(CE_Warning, CPLE_AppDefined,
562 : "Cannot establish CRS because of invalid value "
563 : "for typeOfHorizontalCRS: %d",
564 : nCRSType);
565 0 : proj_destroy(crs);
566 0 : return false;
567 : }
568 :
569 : const char *pszPROJJSON =
570 14 : proj_as_projjson(pjCtxt, crs, nullptr);
571 14 : oSRS.SetFromUserInput(pszPROJJSON);
572 :
573 14 : proj_destroy(crs);
574 : }
575 : else
576 : {
577 0 : CPLError(CE_Warning, CPLE_AppDefined,
578 : "Cannot establish CRS because nameOfHorizontalCRS, "
579 : "typeOfHorizontalCRS, horizontalCS and/or "
580 : "horizontalDatum are missing");
581 0 : return false;
582 : }
583 : }
584 : }
585 : else
586 : {
587 : auto poHorizontalDatumReference =
588 69 : poRootGroup->GetAttribute("horizontalDatumReference");
589 : auto poHorizontalDatumValue =
590 69 : poRootGroup->GetAttribute("horizontalDatumValue");
591 23 : if (poHorizontalDatumReference && poHorizontalDatumValue)
592 : {
593 : const char *pszAuthName =
594 23 : poHorizontalDatumReference->ReadAsString();
595 23 : const char *pszAuthCode = poHorizontalDatumValue->ReadAsString();
596 23 : if (pszAuthName && pszAuthCode)
597 : {
598 46 : if (oSRS.SetFromUserInput(
599 46 : (std::string(pszAuthName) + ':' + pszAuthCode).c_str(),
600 : OGRSpatialReference::
601 23 : SET_FROM_USER_INPUT_LIMITATIONS_get()) !=
602 : OGRERR_NONE)
603 : {
604 0 : oSRS.Clear();
605 : }
606 : }
607 : }
608 : }
609 85 : return !oSRS.IsEmpty();
610 : }
611 :
612 : /************************************************************************/
613 : /* S100GetNumPointsLongitudinalLatitudinal() */
614 : /************************************************************************/
615 :
616 105 : bool S100GetNumPointsLongitudinalLatitudinal(const GDALGroup *poGroup,
617 : int &nNumPointsLongitudinal,
618 : int &nNumPointsLatitudinal)
619 : {
620 315 : auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
621 315 : auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
622 : auto poNumPointsLongitudinal =
623 315 : poGroup->GetAttribute("numPointsLongitudinal");
624 315 : auto poNumPointsLatitudinal = poGroup->GetAttribute("numPointsLatitudinal");
625 105 : if (poSpacingX &&
626 315 : poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
627 210 : poSpacingY &&
628 210 : poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64 &&
629 105 : poNumPointsLongitudinal &&
630 105 : GDALDataTypeIsInteger(
631 210 : poNumPointsLongitudinal->GetDataType().GetNumericDataType()) &&
632 315 : poNumPointsLatitudinal &&
633 105 : GDALDataTypeIsInteger(
634 105 : poNumPointsLatitudinal->GetDataType().GetNumericDataType()))
635 : {
636 105 : nNumPointsLongitudinal = poNumPointsLongitudinal->ReadAsInt();
637 105 : nNumPointsLatitudinal = poNumPointsLatitudinal->ReadAsInt();
638 :
639 : // Those are optional, but use them when available, to detect
640 : // potential inconsistency
641 315 : auto poEastBoundLongitude = poGroup->GetAttribute("eastBoundLongitude");
642 315 : auto poWestBoundLongitude = poGroup->GetAttribute("westBoundLongitude");
643 : auto poSouthBoundLongitude =
644 315 : poGroup->GetAttribute("southBoundLatitude");
645 210 : auto poNorthBoundLatitude = poGroup->GetAttribute("northBoundLatitude");
646 109 : if (poEastBoundLongitude &&
647 4 : GDALDataTypeIsFloating(
648 113 : poEastBoundLongitude->GetDataType().GetNumericDataType()) &&
649 4 : poWestBoundLongitude &&
650 4 : GDALDataTypeIsFloating(
651 8 : poWestBoundLongitude->GetDataType().GetNumericDataType()) &&
652 4 : poSouthBoundLongitude &&
653 4 : GDALDataTypeIsFloating(
654 8 : poSouthBoundLongitude->GetDataType().GetNumericDataType()) &&
655 113 : poNorthBoundLatitude &&
656 4 : GDALDataTypeIsFloating(
657 4 : poNorthBoundLatitude->GetDataType().GetNumericDataType()))
658 : {
659 4 : const double dfSpacingX = poSpacingX->ReadAsDouble();
660 4 : const double dfSpacingY = poSpacingY->ReadAsDouble();
661 :
662 4 : const double dfEast = poEastBoundLongitude->ReadAsDouble();
663 4 : const double dfWest = poWestBoundLongitude->ReadAsDouble();
664 4 : const double dfSouth = poSouthBoundLongitude->ReadAsDouble();
665 4 : const double dfNorth = poNorthBoundLatitude->ReadAsDouble();
666 4 : if (std::fabs((dfWest + nNumPointsLongitudinal * dfSpacingX) -
667 4 : dfEast) < 5 * dfSpacingX &&
668 4 : std::fabs((dfSouth + nNumPointsLatitudinal * dfSpacingY) -
669 4 : dfNorth) < 5 * dfSpacingY)
670 : {
671 : // We need up to 5 spacings for product
672 : // S-111 Trial Data Set Release 1.1/111UK_20210401T000000Z_SolentAndAppr_dcf2.h5
673 : }
674 : else
675 : {
676 0 : CPLError(
677 : CE_Warning, CPLE_AppDefined,
678 : "Caution: "
679 : "eastBoundLongitude/westBoundLongitude/southBoundLatitude/"
680 : "northBoundLatitude/gridSpacingLongitudinal/"
681 : "gridSpacingLatitudinal/numPointsLongitudinal/"
682 : "numPointsLatitudinal do not seem to be consistent");
683 0 : CPLDebug("S100", "Computed east = %f. Actual = %f",
684 0 : dfWest + nNumPointsLongitudinal * dfSpacingX, dfEast);
685 0 : CPLDebug("S100", "Computed north = %f. Actual = %f",
686 0 : dfSouth + nNumPointsLatitudinal * dfSpacingY, dfNorth);
687 : }
688 : }
689 :
690 105 : return true;
691 : }
692 0 : return false;
693 : }
694 :
695 : /************************************************************************/
696 : /* S100GetGeoTransform() */
697 : /************************************************************************/
698 :
699 49 : bool S100GetGeoTransform(const GDALGroup *poGroup, GDALGeoTransform >,
700 : bool bNorthUp)
701 : {
702 147 : auto poOriginX = poGroup->GetAttribute("gridOriginLongitude");
703 147 : auto poOriginY = poGroup->GetAttribute("gridOriginLatitude");
704 147 : auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
705 147 : auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
706 49 : if (poOriginX &&
707 147 : poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 &&
708 98 : poOriginY &&
709 98 : poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 &&
710 98 : poSpacingX &&
711 98 : poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
712 147 : poSpacingY &&
713 49 : poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64)
714 : {
715 49 : int nNumPointsLongitudinal = 0;
716 49 : int nNumPointsLatitudinal = 0;
717 49 : if (!S100GetNumPointsLongitudinalLatitudinal(
718 : poGroup, nNumPointsLongitudinal, nNumPointsLatitudinal))
719 0 : return false;
720 :
721 49 : const double dfSpacingX = poSpacingX->ReadAsDouble();
722 49 : const double dfSpacingY = poSpacingY->ReadAsDouble();
723 :
724 49 : gt[0] = poOriginX->ReadAsDouble();
725 98 : gt[3] = poOriginY->ReadAsDouble() +
726 49 : (bNorthUp ? dfSpacingY * (nNumPointsLatitudinal - 1) : 0);
727 49 : gt[1] = dfSpacingX;
728 49 : gt[5] = bNorthUp ? -dfSpacingY : dfSpacingY;
729 :
730 : // From pixel-center convention to pixel-corner convention
731 49 : gt[0] -= gt[1] / 2;
732 49 : gt[3] -= gt[5] / 2;
733 :
734 49 : return true;
735 : }
736 0 : return false;
737 : }
738 :
739 : /************************************************************************/
740 : /* S100GetDimensions() */
741 : /************************************************************************/
742 :
743 30 : bool S100GetDimensions(
744 : const GDALGroup *poGroup,
745 : std::vector<std::shared_ptr<GDALDimension>> &apoDims,
746 : std::vector<std::shared_ptr<GDALMDArray>> &apoIndexingVars)
747 : {
748 90 : auto poOriginX = poGroup->GetAttribute("gridOriginLongitude");
749 90 : auto poOriginY = poGroup->GetAttribute("gridOriginLatitude");
750 90 : auto poSpacingX = poGroup->GetAttribute("gridSpacingLongitudinal");
751 90 : auto poSpacingY = poGroup->GetAttribute("gridSpacingLatitudinal");
752 30 : if (poOriginX &&
753 90 : poOriginX->GetDataType().GetNumericDataType() == GDT_Float64 &&
754 60 : poOriginY &&
755 60 : poOriginY->GetDataType().GetNumericDataType() == GDT_Float64 &&
756 60 : poSpacingX &&
757 60 : poSpacingX->GetDataType().GetNumericDataType() == GDT_Float64 &&
758 90 : poSpacingY &&
759 30 : poSpacingY->GetDataType().GetNumericDataType() == GDT_Float64)
760 : {
761 30 : int nNumPointsLongitudinal = 0;
762 30 : int nNumPointsLatitudinal = 0;
763 30 : if (!S100GetNumPointsLongitudinalLatitudinal(
764 : poGroup, nNumPointsLongitudinal, nNumPointsLatitudinal))
765 0 : return false;
766 :
767 : {
768 : auto poDim = std::make_shared<GDALDimensionWeakIndexingVar>(
769 60 : std::string(), "Y", GDAL_DIM_TYPE_HORIZONTAL_Y, std::string(),
770 60 : nNumPointsLatitudinal);
771 : auto poIndexingVar = GDALMDArrayRegularlySpaced::Create(
772 90 : std::string(), poDim->GetName(), poDim,
773 120 : poOriginY->ReadAsDouble(), poSpacingY->ReadAsDouble(), 0);
774 30 : poDim->SetIndexingVariable(poIndexingVar);
775 30 : apoDims.emplace_back(poDim);
776 30 : apoIndexingVars.emplace_back(poIndexingVar);
777 : }
778 :
779 : {
780 : auto poDim = std::make_shared<GDALDimensionWeakIndexingVar>(
781 60 : std::string(), "X", GDAL_DIM_TYPE_HORIZONTAL_X, std::string(),
782 60 : nNumPointsLongitudinal);
783 : auto poIndexingVar = GDALMDArrayRegularlySpaced::Create(
784 90 : std::string(), poDim->GetName(), poDim,
785 120 : poOriginX->ReadAsDouble(), poSpacingX->ReadAsDouble(), 0);
786 30 : poDim->SetIndexingVariable(poIndexingVar);
787 30 : apoDims.emplace_back(poDim);
788 30 : apoIndexingVars.emplace_back(poIndexingVar);
789 : }
790 :
791 30 : return true;
792 : }
793 0 : return false;
794 : }
795 :
796 : /************************************************************************/
797 : /* S100ReadVerticalDatum() */
798 : /************************************************************************/
799 :
800 111 : void S100ReadVerticalDatum(GDALMajorObject *poMO, const GDALGroup *poGroup)
801 : {
802 : // https://iho.int/uploads/user/pubs/standards/s-100/S-100_5.2.0_Final_Clean.pdf
803 : // Table 10c-25 - Vertical and sounding datum, page 53
804 : static const struct
805 : {
806 : int nCode;
807 : const char *pszMeaning;
808 : const char *pszAbbrev;
809 : } asVerticalDatums[] = {
810 : {1, "meanLowWaterSprings", "MLWS"},
811 : {2, "meanLowerLowWaterSprings", nullptr},
812 : {3, "meanSeaLevel", "MSL"},
813 : {4, "lowestLowWater", nullptr},
814 : {5, "meanLowWater", "MLW"},
815 : {6, "lowestLowWaterSprings", nullptr},
816 : {7, "approximateMeanLowWaterSprings", nullptr},
817 : {8, "indianSpringLowWater", nullptr},
818 : {9, "lowWaterSprings", nullptr},
819 : {10, "approximateLowestAstronomicalTide", nullptr},
820 : {11, "nearlyLowestLowWater", nullptr},
821 : {12, "meanLowerLowWater", "MLLW"},
822 : {13, "lowWater", "LW"},
823 : {14, "approximateMeanLowWater", nullptr},
824 : {15, "approximateMeanLowerLowWater", nullptr},
825 : {16, "meanHighWater", "MHW"},
826 : {17, "meanHighWaterSprings", "MHWS"},
827 : {18, "highWater", "HW"},
828 : {19, "approximateMeanSeaLevel", nullptr},
829 : {20, "highWaterSprings", nullptr},
830 : {21, "meanHigherHighWater", "MHHW"},
831 : {22, "equinoctialSpringLowWater", nullptr},
832 : {23, "lowestAstronomicalTide", "LAT"},
833 : {24, "localDatum", nullptr},
834 : {25, "internationalGreatLakesDatum1985", nullptr},
835 : {26, "meanWaterLevel", nullptr},
836 : {27, "lowerLowWaterLargeTide", nullptr},
837 : {28, "higherHighWaterLargeTide", nullptr},
838 : {29, "nearlyHighestHighWater", nullptr},
839 : {30, "highestAstronomicalTide", "HAT"},
840 : {44, "balticSeaChartDatum2000", nullptr},
841 : {46, "internationalGreatLakesDatum2020", nullptr},
842 : {47, "seaFloor", nullptr},
843 : {48, "seaSurface", nullptr},
844 : {49, "hydrographicZero", nullptr},
845 : };
846 :
847 111 : int nVerticalDatumReference = 1;
848 : auto poVerticalDatumReference =
849 333 : poGroup->GetAttribute("verticalDatumReference");
850 167 : if (poVerticalDatumReference &&
851 167 : poVerticalDatumReference->GetDataType().GetClass() == GEDTC_NUMERIC)
852 : {
853 56 : nVerticalDatumReference = poVerticalDatumReference->ReadAsInt();
854 : }
855 :
856 333 : auto poVerticalDatum = poGroup->GetAttribute("verticalDatum");
857 178 : if (poVerticalDatum &&
858 178 : poVerticalDatum->GetDataType().GetClass() == GEDTC_NUMERIC)
859 : {
860 67 : if (nVerticalDatumReference == 1)
861 : {
862 67 : bool bFound = false;
863 67 : const auto nVal = poVerticalDatum->ReadAsInt();
864 810 : for (const auto &sVerticalDatum : asVerticalDatums)
865 : {
866 810 : if (sVerticalDatum.nCode == nVal)
867 : {
868 67 : bFound = true;
869 67 : poMO->GDALMajorObject::SetMetadataItem(
870 67 : S100_VERTICAL_DATUM_MEANING, sVerticalDatum.pszMeaning);
871 67 : if (sVerticalDatum.pszAbbrev)
872 67 : poMO->GDALMajorObject::SetMetadataItem(
873 : S100_VERTICAL_DATUM_ABBREV,
874 67 : sVerticalDatum.pszAbbrev);
875 67 : break;
876 : }
877 : }
878 67 : if (!bFound)
879 : {
880 0 : poMO->GDALMajorObject::SetMetadataItem("verticalDatum",
881 : CPLSPrintf("%d", nVal));
882 : }
883 : }
884 0 : else if (nVerticalDatumReference == 2)
885 : {
886 0 : const auto nVal = poVerticalDatum->ReadAsInt();
887 0 : PJ *datum = proj_create_from_database(
888 : OSRGetProjTLSContext(), "EPSG", CPLSPrintf("%d", nVal),
889 : PJ_CATEGORY_DATUM, false, nullptr);
890 0 : poMO->GDALMajorObject::SetMetadataItem("VERTICAL_DATUM_EPSG_CODE",
891 : CPLSPrintf("%d", nVal));
892 0 : if (datum)
893 : {
894 0 : poMO->GDALMajorObject::SetMetadataItem(S100_VERTICAL_DATUM_NAME,
895 : proj_get_name(datum));
896 0 : proj_destroy(datum);
897 : }
898 : }
899 : else
900 : {
901 0 : CPLDebug("S100", "Unhandled verticalDatumReference = %d",
902 : nVerticalDatumReference);
903 : }
904 : }
905 111 : }
906 :
907 : /************************************************************************/
908 : /* S100ReadMetadata() */
909 : /************************************************************************/
910 :
911 55 : std::string S100ReadMetadata(GDALDataset *poDS, const std::string &osFilename,
912 : const GDALGroup *poRootGroup)
913 : {
914 55 : std::string osMetadataFile;
915 747 : for (const auto &poAttr : poRootGroup->GetAttributes())
916 : {
917 692 : const auto &osName = poAttr->GetName();
918 692 : if (osName == "metadata")
919 : {
920 30 : const char *pszVal = poAttr->ReadAsString();
921 30 : if (pszVal && pszVal[0])
922 : {
923 30 : if (CPLHasPathTraversal(pszVal))
924 : {
925 0 : CPLError(CE_Warning, CPLE_AppDefined,
926 : "Path traversal detected in %s", pszVal);
927 0 : continue;
928 : }
929 60 : osMetadataFile = CPLFormFilenameSafe(
930 60 : CPLGetPathSafe(osFilename.c_str()).c_str(), pszVal,
931 30 : nullptr);
932 : VSIStatBufL sStat;
933 30 : if (VSIStatL(osMetadataFile.c_str(), &sStat) != 0)
934 : {
935 : // Test products from https://data.admiralty.co.uk/portal/apps/sites/#/marine-data-portal/pages/s-100
936 : // advertise a metadata filename starting with "MD_", per the spec,
937 : // but the actual filename does not start with "MD_"...
938 6 : if (STARTS_WITH(pszVal, "MD_"))
939 : {
940 12 : osMetadataFile = CPLFormFilenameSafe(
941 12 : CPLGetPathSafe(osFilename.c_str()).c_str(),
942 6 : pszVal + strlen("MD_"), nullptr);
943 6 : if (VSIStatL(osMetadataFile.c_str(), &sStat) != 0)
944 : {
945 6 : osMetadataFile.clear();
946 : }
947 : }
948 : else
949 : {
950 0 : osMetadataFile.clear();
951 : }
952 : }
953 : }
954 : }
955 1280 : else if (osName != "horizontalCRS" &&
956 1225 : osName != "horizontalDatumReference" &&
957 1203 : osName != "horizontalDatumValue" &&
958 1137 : osName != "productSpecification" &&
959 1071 : osName != "eastBoundLongitude" &&
960 1049 : osName != "northBoundLatitude" &&
961 1027 : osName != "southBoundLatitude" &&
962 1502 : osName != "westBoundLongitude" && osName != "extentTypeCode" &&
963 1359 : osName != "verticalCS" && osName != "verticalCoordinateBase" &&
964 774 : osName != "verticalDatumReference" &&
965 971 : osName != "verticalDatum" && osName != "nameOfHorizontalCRS" &&
966 846 : osName != "typeOfHorizontalCRS" && osName != "horizontalCS" &&
967 522 : osName != "horizontalDatum" &&
968 507 : osName != "nameOfHorizontalDatum" &&
969 756 : osName != "primeMeridian" && osName != "spheroid" &&
970 490 : osName != "projectionMethod" &&
971 466 : osName != "projectionParameter1" &&
972 442 : osName != "projectionParameter2" &&
973 421 : osName != "projectionParameter3" &&
974 407 : osName != "projectionParameter4" &&
975 399 : osName != "projectionParameter5" &&
976 1478 : osName != "falseNorthing" && osName != "falseEasting")
977 : {
978 174 : const char *pszVal = poAttr->ReadAsString();
979 174 : if (pszVal && pszVal[0])
980 174 : poDS->GDALDataset::SetMetadataItem(osName.c_str(), pszVal);
981 : }
982 : }
983 55 : return osMetadataFile;
984 : }
|
