Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GRIB Driver
4 : * Purpose: GDALDataset driver for GRIB translator for write support
5 : * Author: Even Rouault <even dot rouault at spatialys dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2017, Even Rouault <even dot rouault at spatialys dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ******************************************************************************
12 : *
13 : */
14 :
15 : /* Support for GRIB2 write capabilities has been funded by Meteorological */
16 : /* Service of Canada */
17 :
18 : #include "cpl_port.h"
19 : #include "gribdataset.h"
20 : #include "gdal_priv_templates.hpp"
21 : #include "memdataset.h"
22 :
23 : #include <algorithm>
24 : #include <cmath>
25 : #include <limits>
26 :
27 : #include "degrib/degrib/meta.h"
28 : CPL_C_START
29 : #include "degrib/g2clib/grib2.h"
30 : CPL_C_END
31 :
32 : /************************************************************************/
33 : /* Lon180to360() */
34 : /************************************************************************/
35 :
36 96 : static inline double Lon180to360(double lon)
37 : {
38 96 : if (lon == 180)
39 0 : return 180;
40 96 : return fmod(fmod(lon, 360) + 360, 360);
41 : }
42 :
43 : /************************************************************************/
44 : /* WriteByte() */
45 : /************************************************************************/
46 :
47 7707 : static bool WriteByte(VSILFILE *fp, int nVal)
48 : {
49 7707 : GByte byVal = static_cast<GByte>(nVal);
50 7707 : return VSIFWriteL(&byVal, 1, sizeof(byVal), fp) == sizeof(byVal);
51 : }
52 :
53 : /************************************************************************/
54 : /* WriteSByte() */
55 : /************************************************************************/
56 :
57 16 : static bool WriteSByte(VSILFILE *fp, int nVal)
58 : {
59 16 : signed char sVal = static_cast<signed char>(nVal);
60 16 : if (sVal == std::numeric_limits<signed char>::min())
61 1 : sVal = std::numeric_limits<signed char>::min() + 1;
62 16 : GByte nUnsignedVal = (sVal < 0) ? static_cast<GByte>(-sVal) | 0x80U
63 : : static_cast<GByte>(sVal);
64 16 : return VSIFWriteL(&nUnsignedVal, 1, sizeof(nUnsignedVal), fp) ==
65 16 : sizeof(nUnsignedVal);
66 : }
67 :
68 : /************************************************************************/
69 : /* WriteUInt16() */
70 : /************************************************************************/
71 :
72 1262 : static bool WriteUInt16(VSILFILE *fp, int nVal)
73 : {
74 1262 : GUInt16 usVal = static_cast<GUInt16>(nVal);
75 1262 : CPL_MSBPTR16(&usVal);
76 1262 : return VSIFWriteL(&usVal, 1, sizeof(usVal), fp) == sizeof(usVal);
77 : }
78 :
79 : /************************************************************************/
80 : /* WriteInt16() */
81 : /************************************************************************/
82 :
83 231 : static bool WriteInt16(VSILFILE *fp, int nVal)
84 : {
85 231 : GInt16 sVal = static_cast<GInt16>(nVal);
86 231 : if (sVal == std::numeric_limits<GInt16>::min())
87 1 : sVal = std::numeric_limits<GInt16>::min() + 1;
88 231 : GUInt16 nUnsignedVal = (sVal < 0) ? static_cast<GUInt16>(-sVal) | 0x8000U
89 : : static_cast<GUInt16>(sVal);
90 231 : CPL_MSBPTR16(&nUnsignedVal);
91 231 : return VSIFWriteL(&nUnsignedVal, 1, sizeof(nUnsignedVal), fp) ==
92 231 : sizeof(nUnsignedVal);
93 : }
94 :
95 : /************************************************************************/
96 : /* WriteUInt32() */
97 : /************************************************************************/
98 :
99 3765 : static bool WriteUInt32(VSILFILE *fp, GUInt32 nVal)
100 : {
101 3765 : CPL_MSBPTR32(&nVal);
102 3765 : return VSIFWriteL(&nVal, 1, sizeof(nVal), fp) == sizeof(nVal);
103 : }
104 :
105 : /************************************************************************/
106 : /* WriteInt32() */
107 : /************************************************************************/
108 :
109 1331 : static bool WriteInt32(VSILFILE *fp, GInt32 nVal)
110 : {
111 1331 : if (nVal == std::numeric_limits<GInt32>::min())
112 0 : nVal = std::numeric_limits<GInt32>::min() + 1;
113 1331 : GUInt32 nUnsignedVal = (nVal < 0)
114 1331 : ? static_cast<GUInt32>(-nVal) | 0x80000000U
115 : : static_cast<GUInt32>(nVal);
116 1331 : CPL_MSBPTR32(&nUnsignedVal);
117 1331 : return VSIFWriteL(&nUnsignedVal, 1, sizeof(nUnsignedVal), fp) ==
118 1331 : sizeof(nUnsignedVal);
119 : }
120 :
121 : /************************************************************************/
122 : /* WriteFloat32() */
123 : /************************************************************************/
124 :
125 200 : static bool WriteFloat32(VSILFILE *fp, float fVal)
126 : {
127 200 : CPL_MSBPTR32(&fVal);
128 200 : return VSIFWriteL(&fVal, 1, sizeof(fVal), fp) == sizeof(fVal);
129 : }
130 :
131 : /************************************************************************/
132 : /* PatchSectionSize() */
133 : /************************************************************************/
134 :
135 348 : static void PatchSectionSize(VSILFILE *fp, vsi_l_offset nStartSection)
136 : {
137 348 : vsi_l_offset nCurOffset = VSIFTellL(fp);
138 348 : VSIFSeekL(fp, nStartSection, SEEK_SET);
139 348 : GUInt32 nSect3Size = static_cast<GUInt32>(nCurOffset - nStartSection);
140 348 : WriteUInt32(fp, nSect3Size);
141 348 : VSIFSeekL(fp, nCurOffset, SEEK_SET);
142 348 : }
143 :
144 : /************************************************************************/
145 : /* GRIB2Section3Writer */
146 : /************************************************************************/
147 :
148 : class GRIB2Section3Writer
149 : {
150 : VSILFILE *fp;
151 : GDALDataset *poSrcDS;
152 : OGRSpatialReference oSRS;
153 : const char *pszProjection;
154 : double dfLLX, dfLLY, dfURX, dfURY;
155 : GDALGeoTransform m_gt{};
156 : int nSplitAndSwapColumn = 0;
157 :
158 : bool WriteScaled(double dfVal, double dfUnit);
159 : bool TransformToGeo(double &dfX, double &dfY);
160 : bool WriteEllipsoidAndRasterSize();
161 :
162 : bool WriteGeographic();
163 : bool WriteRotatedLatLon(double dfLatSouthernPole, double dfLonSouthernPole,
164 : double dfAxisRotation);
165 : bool WriteMercator1SP();
166 : bool WriteMercator2SP(OGRSpatialReference *poSRS = nullptr);
167 : bool WriteTransverseMercator();
168 : bool WritePolarStereographic();
169 : bool WriteLCC1SP();
170 : bool WriteLCC2SPOrAEA(OGRSpatialReference *poSRS = nullptr);
171 : bool WriteLAEA();
172 :
173 : public:
174 : GRIB2Section3Writer(VSILFILE *fpIn, GDALDataset *poSrcDSIn);
175 :
176 162 : inline int SplitAndSwap() const
177 : {
178 162 : return nSplitAndSwapColumn;
179 : }
180 :
181 : bool Write();
182 : };
183 :
184 : /************************************************************************/
185 : /* GRIB2Section3Writer() */
186 : /************************************************************************/
187 :
188 162 : GRIB2Section3Writer::GRIB2Section3Writer(VSILFILE *fpIn, GDALDataset *poSrcDSIn)
189 162 : : fp(fpIn), poSrcDS(poSrcDSIn)
190 : {
191 162 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
192 162 : oSRS.importFromWkt(poSrcDS->GetProjectionRef());
193 162 : pszProjection = oSRS.GetAttrValue("PROJECTION");
194 :
195 162 : poSrcDS->GetGeoTransform(m_gt);
196 :
197 162 : dfLLX = m_gt[0] + m_gt[1] / 2;
198 162 : dfLLY = m_gt[3] + m_gt[5] / 2 + (poSrcDS->GetRasterYSize() - 1) * m_gt[5];
199 162 : dfURX = m_gt[0] + m_gt[1] / 2 + (poSrcDS->GetRasterXSize() - 1) * m_gt[1];
200 162 : dfURY = m_gt[3] + m_gt[5] / 2;
201 162 : if (dfURY < dfLLY)
202 : {
203 0 : double dfTemp = dfURY;
204 0 : dfURY = dfLLY;
205 0 : dfLLY = dfTemp;
206 : }
207 162 : }
208 :
209 : /************************************************************************/
210 : /* WriteEllipsoidAndRasterSize() */
211 : /************************************************************************/
212 :
213 162 : bool GRIB2Section3Writer::WriteEllipsoidAndRasterSize()
214 : {
215 162 : const double dfSemiMajor = oSRS.GetSemiMajor();
216 162 : const double dfSemiMinor = oSRS.GetSemiMinor();
217 162 : const double dfInvFlattening = oSRS.GetInvFlattening();
218 231 : if (std::abs(dfSemiMajor - 6378137.0) < 0.01 &&
219 69 : std::abs(dfInvFlattening - 298.257223563) < 1e-9) // WGS84
220 : {
221 68 : WriteByte(fp, 5); // WGS84
222 68 : WriteByte(fp, GRIB2MISSING_u1);
223 68 : WriteUInt32(fp, GRIB2MISSING_u4);
224 68 : WriteByte(fp, GRIB2MISSING_u1);
225 68 : WriteUInt32(fp, GRIB2MISSING_u4);
226 68 : WriteByte(fp, GRIB2MISSING_u1);
227 68 : WriteUInt32(fp, GRIB2MISSING_u4);
228 : }
229 95 : else if (std::abs(dfSemiMajor - 6378137.0) < 0.01 &&
230 1 : std::abs(dfInvFlattening - 298.257222101) < 1e-9) // GRS80
231 : {
232 1 : WriteByte(fp, 4); // GRS80
233 1 : WriteByte(fp, GRIB2MISSING_u1);
234 1 : WriteUInt32(fp, GRIB2MISSING_u4);
235 1 : WriteByte(fp, GRIB2MISSING_u1);
236 1 : WriteUInt32(fp, GRIB2MISSING_u4);
237 1 : WriteByte(fp, GRIB2MISSING_u1);
238 1 : WriteUInt32(fp, GRIB2MISSING_u4);
239 : }
240 93 : else if (dfInvFlattening == 0)
241 : {
242 : // Earth assumed spherical with radius specified (in m)
243 : // by data producer
244 14 : WriteByte(fp, 1);
245 14 : WriteByte(fp, 2); // scale = * 100
246 14 : WriteUInt32(fp, static_cast<GUInt32>(dfSemiMajor * 100 + 0.5));
247 14 : WriteByte(fp, GRIB2MISSING_u1);
248 14 : WriteUInt32(fp, GRIB2MISSING_u4);
249 14 : WriteByte(fp, GRIB2MISSING_u1);
250 14 : WriteUInt32(fp, GRIB2MISSING_u4);
251 : }
252 : else
253 : {
254 : // Earth assumed oblate spheroid with major and minor axes
255 : // specified (in m) by data producer
256 79 : WriteByte(fp, 7);
257 79 : WriteByte(fp, GRIB2MISSING_u1);
258 79 : WriteUInt32(fp, GRIB2MISSING_u4);
259 79 : WriteByte(fp, 2); // scale = * 100
260 79 : WriteUInt32(fp, static_cast<GUInt32>(dfSemiMajor * 100 + 0.5));
261 79 : WriteByte(fp, 2); // scale = * 100
262 79 : WriteUInt32(fp, static_cast<GUInt32>(dfSemiMinor * 100 + 0.5));
263 : }
264 162 : WriteUInt32(fp, poSrcDS->GetRasterXSize());
265 162 : WriteUInt32(fp, poSrcDS->GetRasterYSize());
266 :
267 162 : return true;
268 : }
269 :
270 : /************************************************************************/
271 : /* WriteScaled() */
272 : /************************************************************************/
273 :
274 1297 : bool GRIB2Section3Writer::WriteScaled(double dfVal, double dfUnit)
275 : {
276 1297 : return WriteInt32(fp, static_cast<GInt32>(floor(dfVal / dfUnit + 0.5)));
277 : }
278 :
279 : /************************************************************************/
280 : /* WriteGeographic() */
281 : /************************************************************************/
282 :
283 71 : bool GRIB2Section3Writer::WriteGeographic()
284 : {
285 71 : WriteUInt16(fp, GS3_LATLON); // Grid template number
286 :
287 71 : WriteEllipsoidAndRasterSize();
288 :
289 77 : if (dfLLX < 0 &&
290 6 : CPLTestBool(CPLGetConfigOption("GRIB_ADJUST_LONGITUDE_RANGE", "YES")))
291 : {
292 6 : CPLDebug("GRIB", "Source longitude range is %lf to %lf", dfLLX, dfURX);
293 6 : double dfOrigLLX = dfLLX;
294 6 : dfLLX = Lon180to360(dfLLX);
295 6 : dfURX = Lon180to360(dfURX);
296 :
297 6 : if (dfLLX > dfURX)
298 : {
299 4 : if (fabs(360 - poSrcDS->GetRasterXSize() * m_gt[1]) < m_gt[1] / 4)
300 : {
301 : // Find the first row number east of the prime meridian
302 4 : nSplitAndSwapColumn =
303 4 : static_cast<int>(ceil((0 - dfOrigLLX) / m_gt[1]));
304 4 : CPLDebug("GRIB",
305 : "Rewrapping around the prime meridian at column %d",
306 : nSplitAndSwapColumn);
307 4 : dfLLX = 0;
308 4 : dfURX = 360 - m_gt[1];
309 : }
310 : else
311 : {
312 0 : CPLDebug("GRIB", "Writing a GRIB with 0-360 longitudes "
313 : "crossing the prime meridian");
314 : }
315 : }
316 6 : CPLDebug("GRIB", "Target longitudes range is %lf %lf", dfLLX, dfURX);
317 : }
318 :
319 71 : WriteUInt32(fp, 0); // Basic angle. 0 equivalent of 1
320 : // Subdivisions of basic angle used. ~0 equivalent of 10^6
321 71 : WriteUInt32(fp, GRIB2MISSING_u4);
322 71 : const double dfAngUnit = 1e-6;
323 71 : WriteScaled(dfLLY, dfAngUnit);
324 71 : WriteScaled(dfLLX, dfAngUnit);
325 71 : WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4); // Resolution and component flags
326 71 : WriteScaled(dfURY, dfAngUnit);
327 71 : WriteScaled(dfURX, dfAngUnit);
328 71 : WriteScaled(m_gt[1], dfAngUnit);
329 71 : WriteScaled(fabs(m_gt[5]), dfAngUnit);
330 71 : WriteByte(fp, GRIB2BIT_2); // Scanning mode: bottom-to-top
331 :
332 71 : return true;
333 : }
334 :
335 : /************************************************************************/
336 : /* WriteRotatedLatLon() */
337 : /************************************************************************/
338 :
339 0 : bool GRIB2Section3Writer::WriteRotatedLatLon(double dfLatSouthernPole,
340 : double dfLonSouthernPole,
341 : double dfAxisRotation)
342 : {
343 0 : WriteUInt16(fp, GS3_ROTATED_LATLON); // Grid template number
344 :
345 0 : WriteEllipsoidAndRasterSize();
346 :
347 0 : if (dfLLX < 0 &&
348 0 : CPLTestBool(CPLGetConfigOption("GRIB_ADJUST_LONGITUDE_RANGE", "YES")))
349 : {
350 0 : CPLDebug("GRIB", "Source longitude range is %lf to %lf", dfLLX, dfURX);
351 0 : double dfOrigLLX = dfLLX;
352 0 : dfLLX = Lon180to360(dfLLX);
353 0 : dfURX = Lon180to360(dfURX);
354 :
355 0 : if (dfLLX > dfURX)
356 : {
357 0 : if (fabs(360 - poSrcDS->GetRasterXSize() * m_gt[1]) < m_gt[1] / 4)
358 : {
359 : // Find the first row number east of the prime meridian
360 0 : nSplitAndSwapColumn =
361 0 : static_cast<int>(ceil((0 - dfOrigLLX) / m_gt[1]));
362 0 : CPLDebug("GRIB",
363 : "Rewrapping around the prime meridian at column %d",
364 : nSplitAndSwapColumn);
365 0 : dfLLX = 0;
366 0 : dfURX = 360 - m_gt[1];
367 : }
368 : else
369 : {
370 0 : CPLDebug("GRIB", "Writing a GRIB with 0-360 longitudes "
371 : "crossing the prime meridian");
372 : }
373 : }
374 0 : CPLDebug("GRIB", "Target longitudes range is %lf %lf", dfLLX, dfURX);
375 : }
376 :
377 0 : WriteUInt32(fp, 0); // Basic angle. 0 equivalent of 1
378 : // Subdivisions of basic angle used. ~0 equivalent of 10^6
379 0 : WriteUInt32(fp, GRIB2MISSING_u4);
380 0 : const double dfAngUnit = 1e-6;
381 0 : WriteScaled(dfLLY, dfAngUnit);
382 0 : WriteScaled(dfLLX, dfAngUnit);
383 0 : WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4); // Resolution and component flags
384 0 : WriteScaled(dfURY, dfAngUnit);
385 0 : WriteScaled(dfURX, dfAngUnit);
386 0 : WriteScaled(m_gt[1], dfAngUnit);
387 0 : WriteScaled(fabs(m_gt[5]), dfAngUnit);
388 0 : WriteByte(fp, GRIB2BIT_2); // Scanning mode: bottom-to-top
389 0 : WriteScaled(dfLatSouthernPole, dfAngUnit);
390 0 : WriteScaled(Lon180to360(dfLonSouthernPole), dfAngUnit);
391 0 : WriteScaled(dfAxisRotation, dfAngUnit);
392 :
393 0 : return true;
394 : }
395 :
396 : /************************************************************************/
397 : /* TransformToGeo() */
398 : /************************************************************************/
399 :
400 21 : bool GRIB2Section3Writer::TransformToGeo(double &dfX, double &dfY)
401 : {
402 42 : OGRSpatialReference oLL; // Construct the "geographic" part of oSRS.
403 21 : oLL.CopyGeogCSFrom(&oSRS);
404 21 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
405 : OGRCoordinateTransformation *poTransformSRSToLL =
406 21 : OGRCreateCoordinateTransformation(&(oSRS), &(oLL));
407 42 : if (poTransformSRSToLL == nullptr ||
408 21 : !poTransformSRSToLL->Transform(1, &dfX, &dfY))
409 : {
410 0 : delete poTransformSRSToLL;
411 0 : return false;
412 : }
413 21 : delete poTransformSRSToLL;
414 21 : if (dfX < 0.0)
415 20 : dfX += 360.0;
416 21 : return true;
417 : }
418 :
419 : /************************************************************************/
420 : /* WriteMercator1SP() */
421 : /************************************************************************/
422 :
423 2 : bool GRIB2Section3Writer::WriteMercator1SP()
424 : {
425 2 : if (oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0) != 0.0)
426 : {
427 0 : CPLError(CE_Failure, CPLE_NotSupported,
428 : "Mercator_1SP with central_meridian != 0 not supported");
429 0 : return false;
430 : }
431 2 : if (oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) != 0.0)
432 : {
433 0 : CPLError(CE_Failure, CPLE_NotSupported,
434 : "Mercator_1SP with latitude_of_origin != 0 not supported");
435 0 : return false;
436 : }
437 :
438 : OGRSpatialReference *poMerc2SP =
439 2 : oSRS.convertToOtherProjection(SRS_PT_MERCATOR_2SP);
440 2 : if (poMerc2SP == nullptr)
441 : {
442 0 : CPLError(CE_Failure, CPLE_NotSupported,
443 : "Cannot get Mercator_2SP formulation");
444 0 : return false;
445 : }
446 :
447 2 : bool bRet = WriteMercator2SP(poMerc2SP);
448 2 : delete poMerc2SP;
449 2 : return bRet;
450 : }
451 :
452 : /************************************************************************/
453 : /* WriteMercator2SP() */
454 : /************************************************************************/
455 :
456 7 : bool GRIB2Section3Writer::WriteMercator2SP(OGRSpatialReference *poSRS)
457 : {
458 7 : if (poSRS == nullptr)
459 5 : poSRS = &oSRS;
460 :
461 7 : if (poSRS->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0) != 0.0)
462 : {
463 0 : CPLError(CE_Failure, CPLE_NotSupported,
464 : "Mercator_2SP with central_meridian != 0 not supported");
465 0 : return false;
466 : }
467 7 : if (poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) != 0.0)
468 : {
469 0 : CPLError(CE_Failure, CPLE_NotSupported,
470 : "Mercator_2SP with latitude_of_origin != 0 not supported");
471 0 : return false;
472 : }
473 :
474 7 : WriteUInt16(fp, GS3_MERCATOR); // Grid template number
475 :
476 7 : WriteEllipsoidAndRasterSize();
477 :
478 7 : if (!TransformToGeo(dfLLX, dfLLY) || !TransformToGeo(dfURX, dfURY))
479 0 : return false;
480 :
481 7 : const double dfAngUnit = 1e-6;
482 7 : WriteScaled(dfLLY, dfAngUnit);
483 7 : WriteScaled(dfLLX, dfAngUnit);
484 7 : WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4); // Resolution and component flags
485 7 : WriteScaled(poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0),
486 : dfAngUnit);
487 7 : WriteScaled(dfURY, dfAngUnit);
488 7 : WriteScaled(dfURX, dfAngUnit);
489 7 : WriteByte(fp, GRIB2BIT_2); // Scanning mode: bottom-to-top
490 7 : WriteInt32(fp, 0); // angle of the grid
491 7 : const double dfLinearUnit = 1e-3;
492 7 : WriteScaled(m_gt[1], dfLinearUnit);
493 7 : WriteScaled(fabs(m_gt[5]), dfLinearUnit);
494 :
495 7 : return true;
496 : }
497 :
498 : /************************************************************************/
499 : /* WriteTransverseMercator() */
500 : /************************************************************************/
501 :
502 78 : bool GRIB2Section3Writer::WriteTransverseMercator()
503 : {
504 78 : WriteUInt16(fp, GS3_TRANSVERSE_MERCATOR); // Grid template number
505 78 : WriteEllipsoidAndRasterSize();
506 :
507 78 : const double dfAngUnit = 1e-6;
508 78 : WriteScaled(oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0),
509 : dfAngUnit);
510 78 : WriteScaled(Lon180to360(oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)),
511 : dfAngUnit);
512 78 : WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4); // Resolution and component flags
513 : float fScale =
514 78 : static_cast<float>(oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0));
515 78 : WriteFloat32(fp, fScale);
516 78 : const double dfLinearUnit = 1e-2;
517 78 : WriteScaled(oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0), dfLinearUnit);
518 78 : WriteScaled(oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0), dfLinearUnit);
519 78 : WriteByte(fp, GRIB2BIT_2); // Scanning mode: bottom-to-top
520 78 : WriteScaled(m_gt[1], dfLinearUnit);
521 78 : WriteScaled(fabs(m_gt[5]), dfLinearUnit);
522 78 : WriteScaled(dfLLX, dfLinearUnit);
523 78 : WriteScaled(dfLLY, dfLinearUnit);
524 78 : WriteScaled(dfURX, dfLinearUnit);
525 78 : WriteScaled(dfURY, dfLinearUnit);
526 :
527 78 : return true;
528 : }
529 :
530 : /************************************************************************/
531 : /* WritePolarStereographic() */
532 : /************************************************************************/
533 :
534 2 : bool GRIB2Section3Writer::WritePolarStereographic()
535 : {
536 2 : WriteUInt16(fp, GS3_POLAR); // Grid template number
537 2 : WriteEllipsoidAndRasterSize();
538 :
539 2 : if (!TransformToGeo(dfLLX, dfLLY))
540 0 : return false;
541 :
542 2 : const double dfAngUnit = 1e-6;
543 2 : WriteScaled(dfLLY, dfAngUnit);
544 2 : WriteScaled(dfLLX, dfAngUnit);
545 2 : WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4); // Resolution and component flags
546 : const double dfLatOrigin =
547 2 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
548 2 : WriteScaled(dfLatOrigin, dfAngUnit);
549 2 : WriteScaled(Lon180to360(oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)),
550 : dfAngUnit);
551 2 : const double dfLinearUnit = 1e-3;
552 2 : WriteScaled(m_gt[1], dfLinearUnit);
553 2 : WriteScaled(fabs(m_gt[5]), dfLinearUnit);
554 : // Projection center flag: BIT1=0 North Pole, BIT1=1 South Pole
555 2 : WriteByte(fp, (dfLatOrigin < 0) ? GRIB2BIT_1 : 0);
556 2 : WriteByte(fp, GRIB2BIT_2); // Scanning mode: bottom-to-top
557 :
558 2 : return true;
559 : }
560 :
561 : /************************************************************************/
562 : /* WriteLCC1SP() */
563 : /************************************************************************/
564 :
565 1 : bool GRIB2Section3Writer::WriteLCC1SP()
566 : {
567 : OGRSpatialReference *poLCC2SP =
568 1 : oSRS.convertToOtherProjection(SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP);
569 1 : if (poLCC2SP == nullptr)
570 : {
571 0 : CPLError(CE_Failure, CPLE_NotSupported,
572 : "Cannot get Lambert_Conformal_Conic_2SP formulation");
573 0 : return false;
574 : }
575 :
576 1 : bool bRet = WriteLCC2SPOrAEA(poLCC2SP);
577 1 : delete poLCC2SP;
578 1 : return bRet;
579 : }
580 :
581 : /************************************************************************/
582 : /* WriteLCC2SPOrAEA() */
583 : /************************************************************************/
584 :
585 3 : bool GRIB2Section3Writer::WriteLCC2SPOrAEA(OGRSpatialReference *poSRS)
586 : {
587 3 : if (poSRS == nullptr)
588 2 : poSRS = &oSRS;
589 3 : if (EQUAL(poSRS->GetAttrValue("PROJECTION"),
590 : SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
591 2 : WriteUInt16(fp, GS3_LAMBERT); // Grid template number
592 : else
593 1 : WriteUInt16(fp, GS3_ALBERS_EQUAL_AREA); // Grid template number
594 :
595 3 : WriteEllipsoidAndRasterSize();
596 :
597 3 : if (!TransformToGeo(dfLLX, dfLLY))
598 0 : return false;
599 :
600 3 : const double dfAngUnit = 1e-6;
601 3 : WriteScaled(dfLLY, dfAngUnit);
602 3 : WriteScaled(dfLLX, dfAngUnit);
603 : // Resolution and component flags. "not applicable" ==> 0 ?
604 3 : WriteByte(fp, 0);
605 3 : WriteScaled(poSRS->GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0),
606 : dfAngUnit);
607 3 : WriteScaled(Lon180to360(oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)),
608 : dfAngUnit);
609 3 : const double dfLinearUnit = 1e-3;
610 3 : WriteScaled(m_gt[1], dfLinearUnit);
611 3 : WriteScaled(fabs(m_gt[5]), dfLinearUnit);
612 3 : WriteByte(fp, 0); // Projection centre flag
613 3 : WriteByte(fp, GRIB2BIT_2); // Scanning mode: bottom-to-top
614 3 : WriteScaled(poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0),
615 : dfAngUnit);
616 3 : WriteScaled(poSRS->GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0),
617 : dfAngUnit);
618 : // Latitude of the southern pole of projection
619 3 : WriteUInt32(fp, GRIB2MISSING_u4);
620 : // Longitude of the southern pole of projection
621 3 : WriteUInt32(fp, GRIB2MISSING_u4);
622 3 : return true;
623 : }
624 :
625 : /************************************************************************/
626 : /* WriteLAEA() */
627 : /************************************************************************/
628 :
629 1 : bool GRIB2Section3Writer::WriteLAEA()
630 : {
631 1 : WriteUInt16(fp, GS3_LAMBERT_AZIMUTHAL); // Grid template number
632 :
633 1 : WriteEllipsoidAndRasterSize();
634 :
635 1 : if (!TransformToGeo(dfLLX, dfLLY) || !TransformToGeo(dfURX, dfURY))
636 0 : return false;
637 :
638 : const bool bNormalizeLongitude =
639 1 : CPLTestBool(CPLGetConfigOption("GRIB_ADJUST_LONGITUDE_RANGE", "YES"));
640 :
641 1 : const double dfAngUnit = 1e-6;
642 1 : WriteScaled(dfLLY, dfAngUnit);
643 1 : if (!bNormalizeLongitude && dfLLX > 360)
644 0 : dfLLX -= 360;
645 1 : WriteScaled(dfLLX, dfAngUnit);
646 1 : WriteScaled(oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0),
647 : dfAngUnit);
648 : const double dfLonCenter =
649 1 : oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
650 1 : WriteScaled(bNormalizeLongitude ? Lon180to360(dfLonCenter) : dfLonCenter,
651 : dfAngUnit);
652 1 : WriteByte(fp, GRIB2BIT_3 | GRIB2BIT_4); // Resolution and component flags
653 1 : const double dfLinearUnit = 1e-3;
654 1 : WriteScaled(m_gt[1], dfLinearUnit);
655 1 : WriteScaled(fabs(m_gt[5]), dfLinearUnit);
656 1 : WriteByte(fp, GRIB2BIT_2); // Scanning mode: bottom-to-top
657 1 : return true;
658 : }
659 :
660 : /************************************************************************/
661 : /* Write() */
662 : /************************************************************************/
663 :
664 162 : bool GRIB2Section3Writer::Write()
665 : {
666 : // Section 3: Grid Definition Section
667 162 : vsi_l_offset nStartSection = VSIFTellL(fp);
668 :
669 162 : WriteUInt32(fp, GRIB2MISSING_u4); // section size
670 :
671 162 : WriteByte(fp, 3); // section number
672 :
673 : // Source of grid definition = Specified in Code Table 3.1
674 162 : WriteByte(fp, 0);
675 :
676 : const GUInt32 nDataPoints =
677 162 : static_cast<GUInt32>(poSrcDS->GetRasterXSize()) *
678 162 : poSrcDS->GetRasterYSize();
679 162 : WriteUInt32(fp, nDataPoints);
680 :
681 : // Number of octets for optional list of numbers defining number of points
682 162 : WriteByte(fp, 0);
683 :
684 : // Interpretation of list of numbers defining number of points =
685 : // No appended list
686 162 : WriteByte(fp, 0);
687 :
688 162 : bool bRet = false;
689 162 : if (oSRS.IsGeographic())
690 : {
691 71 : if (oSRS.IsDerivedGeographic())
692 : {
693 : const OGR_SRSNode *poConversion =
694 0 : oSRS.GetAttrNode("DERIVINGCONVERSION");
695 0 : const char *pszMethod = oSRS.GetAttrValue("METHOD");
696 0 : if (!pszMethod)
697 0 : pszMethod = "unknown";
698 :
699 0 : std::map<std::string, double> oValMap;
700 0 : if (poConversion)
701 : {
702 0 : for (int iChild = 0; iChild < poConversion->GetChildCount();
703 : iChild++)
704 : {
705 0 : const OGR_SRSNode *poNode = poConversion->GetChild(iChild);
706 0 : if (!EQUAL(poNode->GetValue(), "PARAMETER") ||
707 0 : poNode->GetChildCount() <= 2)
708 0 : continue;
709 0 : const char *pszParamStr = poNode->GetChild(0)->GetValue();
710 0 : const char *pszParamVal = poNode->GetChild(1)->GetValue();
711 0 : oValMap[pszParamStr] = CPLAtof(pszParamVal);
712 : }
713 : }
714 :
715 0 : if (poConversion && EQUAL(pszMethod, "PROJ ob_tran o_proj=longlat"))
716 : {
717 0 : const double dfLon0 = oValMap["lon_0"];
718 0 : const double dfLonp = oValMap["o_lon_p"];
719 0 : const double dfLatp = oValMap["o_lat_p"];
720 :
721 0 : const double dfLatSouthernPole = -dfLatp;
722 0 : const double dfLonSouthernPole = dfLon0;
723 0 : const double dfAxisRotation = -dfLonp;
724 0 : bRet = WriteRotatedLatLon(dfLatSouthernPole, dfLonSouthernPole,
725 0 : dfAxisRotation);
726 : }
727 0 : else if (poConversion &&
728 0 : EQUAL(pszMethod, "Pole rotation (netCDF CF convention)"))
729 : {
730 : const double dfGridNorthPoleLat =
731 0 : oValMap["Grid north pole latitude (netCDF CF convention)"];
732 : const double dfGridNorthPoleLong =
733 0 : oValMap["Grid north pole longitude (netCDF CF convention)"];
734 : const double dfNorthPoleGridLong =
735 0 : oValMap["North pole grid longitude (netCDF CF convention)"];
736 :
737 0 : const double dfLon0 = 180.0 + dfGridNorthPoleLong;
738 0 : const double dfLonp = dfNorthPoleGridLong;
739 0 : const double dfLatp = dfGridNorthPoleLat;
740 :
741 0 : const double dfLatSouthernPole = -dfLatp;
742 0 : const double dfLonSouthernPole = dfLon0;
743 0 : const double dfAxisRotation = -dfLonp;
744 0 : bRet = WriteRotatedLatLon(dfLatSouthernPole, dfLonSouthernPole,
745 0 : dfAxisRotation);
746 : }
747 0 : else if (poConversion &&
748 0 : EQUAL(pszMethod, "Pole rotation (GRIB convention)"))
749 : {
750 : const double dfLatSouthernPole =
751 0 : oValMap["Latitude of the southern pole (GRIB convention)"];
752 : const double dfLonSouthernPole =
753 0 : oValMap["Longitude of the southern pole (GRIB convention)"];
754 : const double dfAxisRotation =
755 0 : oValMap["Axis rotation (GRIB convention)"];
756 :
757 0 : bRet = WriteRotatedLatLon(dfLatSouthernPole, dfLonSouthernPole,
758 0 : dfAxisRotation);
759 : }
760 : else
761 : {
762 0 : CPLError(CE_Failure, CPLE_NotSupported,
763 : "Unsupported method for DerivedGeographicCRS: %s",
764 : pszMethod);
765 0 : return false;
766 : }
767 : }
768 : else
769 : {
770 71 : bRet = WriteGeographic();
771 : }
772 : }
773 91 : else if (pszProjection && EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
774 : {
775 2 : bRet = WriteMercator1SP();
776 : }
777 89 : else if (pszProjection && EQUAL(pszProjection, SRS_PT_MERCATOR_2SP))
778 : {
779 5 : bRet = WriteMercator2SP();
780 : }
781 84 : else if (pszProjection && EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
782 : {
783 78 : bRet = WriteTransverseMercator();
784 : }
785 6 : else if (pszProjection && EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
786 : {
787 2 : bRet = WritePolarStereographic();
788 : }
789 4 : else if (pszProjection != nullptr &&
790 4 : EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
791 : {
792 1 : bRet = WriteLCC1SP();
793 : }
794 3 : else if (pszProjection != nullptr &&
795 3 : (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) ||
796 2 : EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA)))
797 : {
798 2 : bRet = WriteLCC2SPOrAEA();
799 : }
800 1 : else if (pszProjection &&
801 1 : EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
802 : {
803 1 : bRet = WriteLAEA();
804 : }
805 :
806 162 : PatchSectionSize(fp, nStartSection);
807 :
808 162 : return bRet;
809 : }
810 :
811 : /************************************************************************/
812 : /* GetBandOption() */
813 : /************************************************************************/
814 :
815 3999 : static const char *GetBandOption(char **papszOptions, GDALDataset *poSrcDS,
816 : int nBand, const char *pszKey,
817 : const char *pszDefault)
818 : {
819 3999 : const char *pszVal = CSLFetchNameValue(
820 : papszOptions, CPLSPrintf("BAND_%d_%s", nBand, pszKey));
821 3999 : if (pszVal == nullptr)
822 : {
823 3998 : pszVal = CSLFetchNameValue(papszOptions, pszKey);
824 : }
825 3999 : if (pszVal == nullptr && poSrcDS != nullptr)
826 : {
827 3134 : pszVal = poSrcDS->GetRasterBand(nBand)->GetMetadataItem(
828 3134 : (CPLString("GRIB_") + pszKey).c_str());
829 : }
830 3999 : if (pszVal == nullptr)
831 : {
832 3613 : pszVal = pszDefault;
833 : }
834 3999 : return pszVal;
835 : }
836 :
837 : /************************************************************************/
838 : /* GRIB2Section567Writer */
839 : /************************************************************************/
840 :
841 : class GRIB2Section567Writer
842 : {
843 : VSILFILE *m_fp;
844 : GDALDataset *m_poSrcDS;
845 : int m_nBand;
846 : int m_nXSize;
847 : int m_nYSize;
848 : GUInt32 m_nDataPoints;
849 : GDALDataType m_eDT;
850 : GDALGeoTransform m_gt{};
851 : int m_nDecimalScaleFactor;
852 : double m_dfDecimalScale;
853 : float m_fMin;
854 : float m_fMax;
855 : double m_dfMinScaled;
856 : int m_nBits;
857 : bool m_bUseZeroBits;
858 : float m_fValOffset;
859 : int m_bHasNoData;
860 : double m_dfNoData;
861 : int m_nSplitAndSwap;
862 :
863 : float *GetFloatData();
864 : bool WriteSimplePacking();
865 : bool WriteComplexPacking(int nSpatialDifferencingOrder);
866 : bool WriteIEEE(GDALProgressFunc pfnProgress, void *pProgressData);
867 : bool WritePNG();
868 : bool WriteJPEG2000(char **papszOptions);
869 :
870 : public:
871 : GRIB2Section567Writer(VSILFILE *fp, GDALDataset *poSrcDS, int nBand,
872 : int nSplitAndSwap);
873 :
874 : bool Write(float fValOffset, char **papszOptions,
875 : GDALProgressFunc pfnProgress, void *pProgressData);
876 : void WriteComplexPackingNoData();
877 : };
878 :
879 : /************************************************************************/
880 : /* GRIB2Section567Writer() */
881 : /************************************************************************/
882 :
883 156 : GRIB2Section567Writer::GRIB2Section567Writer(VSILFILE *fp, GDALDataset *poSrcDS,
884 156 : int nBand, int nSplitAndSwap)
885 : : m_fp(fp), m_poSrcDS(poSrcDS), m_nBand(nBand),
886 468 : m_nXSize(poSrcDS->GetRasterXSize()), m_nYSize(poSrcDS->GetRasterYSize()),
887 156 : m_nDataPoints(static_cast<GUInt32>(m_nXSize) * m_nYSize),
888 156 : m_eDT(m_poSrcDS->GetRasterBand(m_nBand)->GetRasterDataType()),
889 : m_nDecimalScaleFactor(0), m_dfDecimalScale(1.0), m_fMin(0.0f),
890 : m_fMax(0.0f), m_dfMinScaled(0.0), m_nBits(0), m_bUseZeroBits(false),
891 : m_fValOffset(0.0), m_bHasNoData(false), m_dfNoData(0.0),
892 312 : m_nSplitAndSwap(nSplitAndSwap)
893 : {
894 156 : m_poSrcDS->GetGeoTransform(m_gt);
895 156 : m_dfNoData = m_poSrcDS->GetRasterBand(nBand)->GetNoDataValue(&m_bHasNoData);
896 156 : }
897 :
898 : /************************************************************************/
899 : /* GetFloatData() */
900 : /************************************************************************/
901 :
902 118 : float *GRIB2Section567Writer::GetFloatData()
903 : {
904 : float *pafData =
905 118 : static_cast<float *>(VSI_MALLOC2_VERBOSE(m_nDataPoints, sizeof(float)));
906 118 : if (pafData == nullptr)
907 : {
908 0 : return nullptr;
909 : }
910 118 : CPLErr eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
911 118 : GF_Read, m_nSplitAndSwap, 0, m_nXSize - m_nSplitAndSwap, m_nYSize,
912 118 : pafData + (m_gt[5] < 0 ? (m_nYSize - 1) * m_nXSize : 0),
913 118 : m_nXSize - m_nSplitAndSwap, m_nYSize, GDT_Float32, sizeof(float),
914 118 : m_gt[5] < 0 ? -static_cast<GSpacing>(m_nXSize * sizeof(float))
915 0 : : static_cast<GSpacing>(m_nXSize * sizeof(float)),
916 : nullptr);
917 118 : if (eErr != CE_None)
918 : {
919 0 : VSIFree(pafData);
920 0 : return nullptr;
921 : }
922 118 : if (m_nSplitAndSwap > 0)
923 : {
924 1 : eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
925 : GF_Read, 0, 0, m_nSplitAndSwap, m_nYSize,
926 1 : pafData + (m_gt[5] < 0 ? (m_nYSize - 1) * m_nXSize : 0) +
927 1 : (m_nXSize - m_nSplitAndSwap),
928 : m_nSplitAndSwap, m_nYSize, GDT_Float32, sizeof(float),
929 1 : m_gt[5] < 0 ? -static_cast<GSpacing>(m_nXSize * sizeof(float))
930 0 : : static_cast<GSpacing>(m_nXSize * sizeof(float)),
931 : nullptr);
932 1 : if (eErr != CE_None)
933 : {
934 0 : VSIFree(pafData);
935 0 : return nullptr;
936 : }
937 : }
938 :
939 118 : m_fMin = std::numeric_limits<float>::max();
940 118 : m_fMax = -std::numeric_limits<float>::max();
941 118 : bool bHasNoDataValuePoint = false;
942 118 : bool bHasDataValuePoint = false;
943 1681620 : for (GUInt32 i = 0; i < m_nDataPoints; i++)
944 : {
945 1681500 : if (m_bHasNoData && pafData[i] == static_cast<float>(m_dfNoData))
946 : {
947 680030 : if (!bHasNoDataValuePoint)
948 8 : bHasNoDataValuePoint = true;
949 680030 : continue;
950 : }
951 1001470 : if (!std::isfinite(pafData[i]))
952 : {
953 0 : CPLError(CE_Failure, CPLE_NotSupported,
954 : "Non-finite values not supported for "
955 : "this data encoding");
956 0 : VSIFree(pafData);
957 0 : return nullptr;
958 : }
959 1001470 : if (!bHasDataValuePoint)
960 116 : bHasDataValuePoint = true;
961 1001470 : pafData[i] += m_fValOffset;
962 1001470 : if (pafData[i] < m_fMin)
963 623 : m_fMin = pafData[i];
964 1001470 : if (pafData[i] > m_fMax)
965 290 : m_fMax = pafData[i];
966 : }
967 118 : if (m_fMin > m_fMax)
968 : {
969 2 : m_fMin = m_fMax = static_cast<float>(m_dfNoData);
970 : }
971 :
972 : // We check that the actual range of values got from the above RasterIO
973 : // request does not go over the expected range of the datatype, as we
974 : // later assume that for computing nMaxBitsPerElt.
975 : // This shouldn't happen for well-behaved drivers, but this can still
976 : // happen in practice, if some drivers don't completely fill buffers etc.
977 177 : if (m_fMax > m_fMin && GDALDataTypeIsInteger(m_eDT) &&
978 59 : ceil(log(m_fMax - m_fMin) / log(2.0)) > GDALGetDataTypeSizeBits(m_eDT))
979 : {
980 0 : CPLError(CE_Failure, CPLE_AppDefined,
981 : "Garbage values found when requesting input dataset");
982 0 : VSIFree(pafData);
983 0 : return nullptr;
984 : }
985 :
986 118 : m_dfMinScaled =
987 118 : m_dfDecimalScale == 1.0 ? m_fMin : floor(m_fMin * m_dfDecimalScale);
988 236 : if (!(m_dfMinScaled >= -std::numeric_limits<float>::max() &&
989 118 : m_dfMinScaled < std::numeric_limits<float>::max()))
990 : {
991 0 : CPLError(CE_Failure, CPLE_AppDefined,
992 : "Scaled min value not representable on IEEE754 "
993 : "single precision float");
994 0 : VSIFree(pafData);
995 0 : return nullptr;
996 : }
997 :
998 118 : const double dfScaledMaxDiff = (m_fMax - m_fMin) * m_dfDecimalScale;
999 118 : if (GDALDataTypeIsFloating(m_eDT) && m_nBits == 0 && dfScaledMaxDiff > 0 &&
1000 : dfScaledMaxDiff <= 256)
1001 : {
1002 8 : m_nBits = 8;
1003 : }
1004 :
1005 118 : m_bUseZeroBits =
1006 194 : (m_fMin == m_fMax && !(bHasDataValuePoint && bHasNoDataValuePoint)) ||
1007 76 : (!GDALDataTypeIsFloating(m_eDT) && dfScaledMaxDiff < 1.0);
1008 :
1009 118 : return pafData;
1010 : }
1011 :
1012 : /************************************************************************/
1013 : /* WriteSimplePacking() */
1014 : /************************************************************************/
1015 :
1016 : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-0.shtml
1017 66 : bool GRIB2Section567Writer::WriteSimplePacking()
1018 : {
1019 66 : float *pafData = GetFloatData();
1020 66 : if (pafData == nullptr)
1021 0 : return false;
1022 :
1023 : const int nBitCorrectionForDec =
1024 66 : static_cast<int>(ceil(m_nDecimalScaleFactor * log(10.0) / log(2.0)));
1025 : const int nMaxBitsPerElt = std::max(
1026 66 : 1, std::min(31, (m_nBits > 0) ? m_nBits
1027 56 : : GDALGetDataTypeSizeBits(m_eDT) +
1028 66 : nBitCorrectionForDec));
1029 66 : if (nMaxBitsPerElt > 0 &&
1030 66 : m_nDataPoints > static_cast<GUInt32>(INT_MAX) / nMaxBitsPerElt)
1031 : {
1032 0 : CPLError(CE_Failure, CPLE_AppDefined,
1033 : "Int overflow while computing maximum number of bits");
1034 0 : VSIFree(pafData);
1035 0 : return false;
1036 : }
1037 :
1038 66 : const int nMaxSize = (m_nDataPoints * nMaxBitsPerElt + 7) / 8;
1039 66 : void *pabyData = VSI_MALLOC_VERBOSE(nMaxSize);
1040 66 : if (pabyData == nullptr)
1041 : {
1042 0 : VSIFree(pafData);
1043 0 : VSIFree(pabyData);
1044 0 : return false;
1045 : }
1046 :
1047 : // Indices expected by simpack()
1048 : enum
1049 : {
1050 : TMPL5_R_IDX = 0, // Reference value (R)
1051 : TMPL5_E_IDX = 1, // Binary scale factor (E)
1052 : TMPL5_D_IDX = 2, // Decimal scale factor (D)
1053 : TMPL5_NBITS_IDX = 3, // Number of bits used for each packed value
1054 : TMPL5_TYPE_IDX = 4 // type of original data
1055 : };
1056 :
1057 66 : g2int idrstmpl[TMPL5_TYPE_IDX + 1] = {0};
1058 66 : idrstmpl[TMPL5_R_IDX] = 0; // reference value, to be filled by simpack
1059 66 : idrstmpl[TMPL5_E_IDX] = 0; // binary scale factor, to be filled by simpack
1060 66 : idrstmpl[TMPL5_D_IDX] = m_nDecimalScaleFactor;
1061 : // to be filled by simpack if set to 0
1062 66 : idrstmpl[TMPL5_NBITS_IDX] = m_nBits;
1063 : // to be filled by simpack (and we will ignore it)
1064 66 : idrstmpl[TMPL5_TYPE_IDX] = 0;
1065 66 : g2int nLengthPacked = 0;
1066 66 : simpack(pafData, m_nDataPoints, idrstmpl,
1067 : static_cast<unsigned char *>(pabyData), &nLengthPacked);
1068 66 : CPLAssert(nLengthPacked <= nMaxSize);
1069 66 : if (nLengthPacked < 0)
1070 : {
1071 0 : CPLError(CE_Failure, CPLE_AppDefined, "Error while packing");
1072 0 : VSIFree(pafData);
1073 0 : VSIFree(pabyData);
1074 0 : return false;
1075 : }
1076 :
1077 : // Section 5: Data Representation Section
1078 66 : WriteUInt32(m_fp, 21); // section size
1079 66 : WriteByte(m_fp, 5); // section number
1080 66 : WriteUInt32(m_fp, m_nDataPoints);
1081 66 : WriteUInt16(m_fp, GS5_SIMPLE);
1082 : float fRefValue;
1083 66 : memcpy(&fRefValue, &idrstmpl[TMPL5_R_IDX], 4);
1084 66 : WriteFloat32(m_fp, fRefValue);
1085 66 : WriteInt16(m_fp, idrstmpl[TMPL5_E_IDX]);
1086 66 : WriteInt16(m_fp, idrstmpl[TMPL5_D_IDX]);
1087 66 : WriteByte(m_fp, idrstmpl[TMPL5_NBITS_IDX]);
1088 : // Type of original data: 0=Floating, 1=Integer
1089 66 : WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
1090 :
1091 : // Section 6: Bitmap section
1092 : #ifdef DEBUG
1093 66 : if (CPLTestBool(CPLGetConfigOption("GRIB_WRITE_BITMAP_TEST", "NO")))
1094 : {
1095 : // Just for the purpose of generating a test product !
1096 : static int counter = 0;
1097 0 : counter++;
1098 0 : if (counter == 1)
1099 : {
1100 0 : WriteUInt32(m_fp, 6 + ((m_nDataPoints + 7) / 8)); // section size
1101 0 : WriteByte(m_fp, 6); // section number
1102 0 : WriteByte(m_fp, 0); // bitmap
1103 0 : for (GUInt32 i = 0; i < (m_nDataPoints + 7) / 8; i++)
1104 0 : WriteByte(m_fp, 255);
1105 : }
1106 : else
1107 : {
1108 0 : WriteUInt32(m_fp, 6); // section size
1109 0 : WriteByte(m_fp, 6); // section number
1110 0 : WriteByte(m_fp, 254); // reuse previous bitmap
1111 : }
1112 : }
1113 : else
1114 : #endif
1115 : {
1116 66 : WriteUInt32(m_fp, 6); // section size
1117 66 : WriteByte(m_fp, 6); // section number
1118 66 : WriteByte(m_fp, GRIB2MISSING_u1); // no bitmap
1119 : }
1120 :
1121 : // Section 7: Data Section
1122 66 : WriteUInt32(m_fp, 5 + nLengthPacked); // section size
1123 66 : WriteByte(m_fp, 7); // section number
1124 66 : if (static_cast<int>(VSIFWriteL(pabyData, 1, nLengthPacked, m_fp)) !=
1125 : nLengthPacked)
1126 : {
1127 0 : VSIFree(pafData);
1128 0 : VSIFree(pabyData);
1129 0 : return false;
1130 : }
1131 :
1132 66 : VSIFree(pafData);
1133 66 : VSIFree(pabyData);
1134 :
1135 66 : return true;
1136 : }
1137 :
1138 : /************************************************************************/
1139 : /* WriteComplexPackingNoData() */
1140 : /************************************************************************/
1141 :
1142 20 : void GRIB2Section567Writer::WriteComplexPackingNoData()
1143 : {
1144 20 : if (!m_bHasNoData)
1145 : {
1146 9 : WriteUInt32(m_fp, GRIB2MISSING_u4);
1147 : }
1148 11 : else if (GDALDataTypeIsFloating(m_eDT))
1149 : {
1150 7 : WriteFloat32(m_fp, static_cast<float>(m_dfNoData));
1151 : }
1152 : else
1153 : {
1154 4 : if (GDALIsValueInRange<int>(m_dfNoData))
1155 : {
1156 4 : WriteInt32(m_fp, static_cast<int>(m_dfNoData));
1157 : }
1158 : else
1159 : {
1160 0 : WriteUInt32(m_fp, GRIB2MISSING_u4);
1161 : }
1162 : }
1163 20 : }
1164 :
1165 : /************************************************************************/
1166 : /* WriteComplexPacking() */
1167 : /************************************************************************/
1168 :
1169 : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-2.shtml
1170 : // and http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-3.shtml
1171 :
1172 22 : bool GRIB2Section567Writer::WriteComplexPacking(int nSpatialDifferencingOrder)
1173 : {
1174 22 : if (nSpatialDifferencingOrder < 0 || nSpatialDifferencingOrder > 2)
1175 : {
1176 1 : CPLError(CE_Failure, CPLE_AppDefined,
1177 : "Unsupported value for SPATIAL_DIFFERENCING_ORDER");
1178 1 : return false;
1179 : }
1180 :
1181 21 : float *pafData = GetFloatData();
1182 21 : if (pafData == nullptr)
1183 0 : return false;
1184 :
1185 21 : const float fNoData = static_cast<float>(m_dfNoData);
1186 21 : if (m_bUseZeroBits)
1187 : {
1188 : // Case where all values are at nodata or a single value
1189 8 : VSIFree(pafData);
1190 :
1191 : // Section 5: Data Representation Section
1192 8 : WriteUInt32(m_fp, 47); // section size
1193 8 : WriteByte(m_fp, 5); // section number
1194 8 : WriteUInt32(m_fp, m_nDataPoints);
1195 8 : WriteUInt16(m_fp, GS5_CMPLX);
1196 8 : WriteFloat32(m_fp, m_fMin); // ref value = nodata or single data
1197 8 : WriteInt16(m_fp, 0); // binary scale factor
1198 8 : WriteInt16(m_fp, 0); // decimal scale factor
1199 8 : WriteByte(m_fp, 0); // number of bits
1200 : // Type of original data: 0=Floating, 1=Integer
1201 8 : WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
1202 8 : WriteByte(m_fp, 0);
1203 8 : WriteByte(m_fp, m_bHasNoData ? 1 : 0); // 1 missing value
1204 8 : WriteComplexPackingNoData();
1205 8 : WriteUInt32(m_fp, GRIB2MISSING_u4);
1206 8 : WriteUInt32(m_fp, 0);
1207 8 : WriteByte(m_fp, 0);
1208 8 : WriteByte(m_fp, 0);
1209 8 : WriteUInt32(m_fp, 0);
1210 8 : WriteByte(m_fp, 0);
1211 8 : WriteUInt32(m_fp, 0);
1212 8 : WriteByte(m_fp, 0);
1213 :
1214 : // Section 6: Bitmap section
1215 8 : WriteUInt32(m_fp, 6); // section size
1216 8 : WriteByte(m_fp, 6); // section number
1217 8 : WriteByte(m_fp, GRIB2MISSING_u1); // no bitmap
1218 :
1219 : // Section 7: Data Section
1220 8 : WriteUInt32(m_fp, 5); // section size
1221 8 : WriteByte(m_fp, 7); // section number
1222 :
1223 8 : return true;
1224 : }
1225 :
1226 : const int nBitCorrectionForDec =
1227 13 : static_cast<int>(ceil(m_nDecimalScaleFactor * log(10.0) / log(2.0)));
1228 : const int nMaxBitsPerElt = std::max(
1229 13 : 1, std::min(31, (m_nBits > 0) ? m_nBits
1230 6 : : GDALGetDataTypeSizeBits(m_eDT) +
1231 13 : nBitCorrectionForDec));
1232 13 : if (nMaxBitsPerElt > 0 &&
1233 13 : m_nDataPoints > static_cast<GUInt32>(INT_MAX) / nMaxBitsPerElt)
1234 : {
1235 0 : CPLError(CE_Failure, CPLE_AppDefined,
1236 : "Int overflow while computing maximum number of bits");
1237 0 : VSIFree(pafData);
1238 0 : return false;
1239 : }
1240 :
1241 : // No idea what is the exact maximum bound... Take the value of simple
1242 : // packing and multiply by 2, plus some constant...
1243 13 : const int nMaxSize = 10000 + 2 * ((m_nDataPoints * nMaxBitsPerElt + 7) / 8);
1244 13 : void *pabyData = VSI_MALLOC_VERBOSE(nMaxSize);
1245 13 : if (pabyData == nullptr)
1246 : {
1247 0 : VSIFree(pafData);
1248 0 : VSIFree(pabyData);
1249 0 : return false;
1250 : }
1251 :
1252 13 : const double dfScaledMaxDiff =
1253 13 : (m_fMax == m_fMin) ? 1 : (m_fMax - m_fMin) * m_dfDecimalScale;
1254 13 : if (m_nBits == 0)
1255 : {
1256 6 : double dfTemp = log(ceil(dfScaledMaxDiff)) / log(2.0);
1257 6 : m_nBits = std::max(1, std::min(31, static_cast<int>(ceil(dfTemp))));
1258 : }
1259 13 : const int nMaxNum = (m_nBits == 31) ? INT_MAX : ((1 << m_nBits) - 1);
1260 13 : double dfTemp = log(nMaxNum / dfScaledMaxDiff) / log(2.0);
1261 13 : int nBinaryScaleFactor = static_cast<GInt16>(ceil(-dfTemp));
1262 :
1263 : // Indices expected by cmplxpack()
1264 : enum
1265 : {
1266 : TMPL5_R_IDX = 0, // reference value
1267 : TMPL5_E_IDX = 1, // binary scale factor
1268 : TMPL5_D_IDX = 2, // decimal scale factor
1269 : TMPL5_NBITS_IDX = 3, // number of bits
1270 : TMPL5_TYPE_IDX = 4, // type of original data
1271 : TMPL5_GROUP_SPLITTING_IDX = 5, // Group splitting method used
1272 : TMPL5_MISSING_VALUE_MGNT_IDX = 6, // Missing value management used
1273 : TMPL5_PRIMARY_MISSING_VALUE_IDX = 7, // Primary missing value
1274 : TMPL5_SECONDARY_MISSING_VALUE_IDX = 8, // Secondary missing value
1275 : TMPL5_NG_IDX = 9, // number of groups of data values
1276 : TMPL5_REF_GROUP_WIDTHS_IDX = 10, // Reference for group widths
1277 : // Number of bits used for the group widths
1278 : TMPL5_NBITS_GROUP_WIDTHS_IDX = 11,
1279 : TMPL5_REF_GROUP_LENGTHS_IDX = 12, // Reference for group lengths
1280 : // Length increment for the group lengths
1281 : TMPL5_LENGTH_INCR_GROUP_LENGTHS_IDX = 13,
1282 : TMPL5_TRUE_LENGTH_LAST_GROUP_IDX = 14, // True length of last group
1283 : // Number of bits used for the scaled group lengths
1284 : TMPL5_NBITS_SCALED_GROUP_LENGTHS_IDX = 15,
1285 : // Order of spatial differencing
1286 : TMPL5_ORDER_SPATIAL_DIFFERENCE_IDX = 16,
1287 : // Number of octets required in the data section to specify extra
1288 : // descriptors needed for spatial differencing
1289 : TMPL5_NB_OCTETS_EXTRA_DESCR_IDX = 17
1290 : };
1291 :
1292 13 : g2int idrstmpl[TMPL5_NB_OCTETS_EXTRA_DESCR_IDX + 1] = {0};
1293 13 : idrstmpl[TMPL5_E_IDX] = nBinaryScaleFactor;
1294 13 : idrstmpl[TMPL5_D_IDX] = m_nDecimalScaleFactor;
1295 13 : idrstmpl[TMPL5_MISSING_VALUE_MGNT_IDX] = m_bHasNoData ? 1 : 0;
1296 13 : idrstmpl[TMPL5_ORDER_SPATIAL_DIFFERENCE_IDX] = nSpatialDifferencingOrder;
1297 13 : if (m_bHasNoData)
1298 : {
1299 4 : memcpy(&idrstmpl[TMPL5_PRIMARY_MISSING_VALUE_IDX], &fNoData, 4);
1300 : }
1301 13 : g2int nLengthPacked = 0;
1302 13 : const int nTemplateNumber =
1303 13 : (nSpatialDifferencingOrder > 0) ? GS5_CMPLXSEC : GS5_CMPLX;
1304 13 : cmplxpack(pafData, m_nDataPoints, nTemplateNumber, idrstmpl,
1305 : static_cast<unsigned char *>(pabyData), &nLengthPacked);
1306 13 : CPLAssert(nLengthPacked <= nMaxSize);
1307 13 : if (nLengthPacked < 0)
1308 : {
1309 1 : CPLError(CE_Failure, CPLE_AppDefined, "Error while packing");
1310 1 : VSIFree(pafData);
1311 1 : VSIFree(pabyData);
1312 1 : return false;
1313 : }
1314 :
1315 : // Section 5: Data Representation Section
1316 12 : WriteUInt32(m_fp, nTemplateNumber == GS5_CMPLX ? 47 : 49); // section size
1317 12 : WriteByte(m_fp, 5); // section number
1318 12 : WriteUInt32(m_fp, m_nDataPoints);
1319 12 : WriteUInt16(m_fp, nTemplateNumber);
1320 : float fRefValue;
1321 12 : memcpy(&fRefValue, &idrstmpl[TMPL5_R_IDX], 4);
1322 12 : WriteFloat32(m_fp, fRefValue);
1323 12 : WriteInt16(m_fp, idrstmpl[TMPL5_E_IDX]);
1324 12 : WriteInt16(m_fp, idrstmpl[TMPL5_D_IDX]);
1325 12 : WriteByte(m_fp, idrstmpl[TMPL5_NBITS_IDX]);
1326 : // Type of original data: 0=Floating, 1=Integer
1327 12 : WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
1328 12 : WriteByte(m_fp, idrstmpl[TMPL5_GROUP_SPLITTING_IDX]);
1329 12 : WriteByte(m_fp, idrstmpl[TMPL5_MISSING_VALUE_MGNT_IDX]);
1330 12 : WriteComplexPackingNoData();
1331 12 : WriteUInt32(m_fp, GRIB2MISSING_u4);
1332 12 : WriteUInt32(m_fp, idrstmpl[TMPL5_NG_IDX]);
1333 12 : WriteByte(m_fp, idrstmpl[TMPL5_REF_GROUP_WIDTHS_IDX]);
1334 12 : WriteByte(m_fp, idrstmpl[TMPL5_NBITS_GROUP_WIDTHS_IDX]);
1335 12 : WriteUInt32(m_fp, idrstmpl[TMPL5_REF_GROUP_LENGTHS_IDX]);
1336 12 : WriteByte(m_fp, idrstmpl[TMPL5_LENGTH_INCR_GROUP_LENGTHS_IDX]);
1337 12 : WriteUInt32(m_fp, idrstmpl[TMPL5_TRUE_LENGTH_LAST_GROUP_IDX]);
1338 12 : WriteByte(m_fp, idrstmpl[TMPL5_NBITS_SCALED_GROUP_LENGTHS_IDX]);
1339 12 : if (nTemplateNumber == GS5_CMPLXSEC)
1340 : {
1341 3 : WriteByte(m_fp, nSpatialDifferencingOrder);
1342 3 : WriteByte(m_fp, idrstmpl[TMPL5_NB_OCTETS_EXTRA_DESCR_IDX]);
1343 : }
1344 :
1345 : // Section 6: Bitmap section
1346 12 : WriteUInt32(m_fp, 6); // section size
1347 12 : WriteByte(m_fp, 6); // section number
1348 12 : WriteByte(m_fp, GRIB2MISSING_u1); // no bitmap
1349 :
1350 : // Section 7: Data Section
1351 12 : WriteUInt32(m_fp, 5 + nLengthPacked); // section size
1352 12 : WriteByte(m_fp, 7); // section number
1353 12 : if (static_cast<int>(VSIFWriteL(pabyData, 1, nLengthPacked, m_fp)) !=
1354 : nLengthPacked)
1355 : {
1356 0 : VSIFree(pafData);
1357 0 : VSIFree(pabyData);
1358 0 : return false;
1359 : }
1360 :
1361 12 : VSIFree(pafData);
1362 12 : VSIFree(pabyData);
1363 :
1364 12 : return true;
1365 : }
1366 :
1367 : /************************************************************************/
1368 : /* WriteIEEE() */
1369 : /************************************************************************/
1370 :
1371 : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-4.shtml
1372 30 : bool GRIB2Section567Writer::WriteIEEE(GDALProgressFunc pfnProgress,
1373 : void *pProgressData)
1374 : {
1375 : GDALDataType eReqDT;
1376 30 : if (GDALGetDataTypeSizeBytes(m_eDT) <= 2 || m_eDT == GDT_Float32)
1377 11 : eReqDT = GDT_Float32;
1378 : else
1379 19 : eReqDT = GDT_Float64;
1380 :
1381 : // Section 5: Data Representation Section
1382 30 : WriteUInt32(m_fp, 12); // section size
1383 30 : WriteByte(m_fp, 5); // section number
1384 30 : WriteUInt32(m_fp, m_nDataPoints);
1385 30 : WriteUInt16(m_fp, GS5_IEEE);
1386 30 : WriteByte(m_fp, (eReqDT == GDT_Float32) ? 1 : 2); // Precision
1387 :
1388 : // Section 6: Bitmap section
1389 30 : WriteUInt32(m_fp, 6); // section size
1390 30 : WriteByte(m_fp, 6); // section number
1391 30 : WriteByte(m_fp, GRIB2MISSING_u1); // no bitmap
1392 :
1393 : // Section 7: Data Section
1394 : const size_t nBufferSize =
1395 30 : static_cast<size_t>(m_nXSize) * GDALGetDataTypeSizeBytes(eReqDT);
1396 : // section size
1397 30 : WriteUInt32(m_fp, static_cast<GUInt32>(5 + nBufferSize * m_nYSize));
1398 30 : WriteByte(m_fp, 7); // section number
1399 30 : void *pData = CPLMalloc(nBufferSize);
1400 30 : const int nDenominator = std::max(1, m_poSrcDS->GetRasterCount());
1401 60 : void *pScaledProgressData = GDALCreateScaledProgress(
1402 30 : static_cast<double>(m_nBand - 1) / nDenominator,
1403 30 : static_cast<double>(m_nBand) / nDenominator, pfnProgress,
1404 : pProgressData);
1405 486 : for (int i = 0; i < m_nYSize; i++)
1406 : {
1407 456 : int iSrcLine = m_gt[5] < 0 ? m_nYSize - 1 - i : i;
1408 456 : CPLErr eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
1409 456 : GF_Read, m_nSplitAndSwap, iSrcLine, m_nXSize - m_nSplitAndSwap, 1,
1410 456 : pData, m_nXSize - m_nSplitAndSwap, 1, eReqDT, 0, 0, nullptr);
1411 456 : if (eErr != CE_None)
1412 : {
1413 0 : CPLFree(pData);
1414 0 : GDALDestroyScaledProgress(pScaledProgressData);
1415 0 : return false;
1416 : }
1417 456 : if (m_nSplitAndSwap > 0)
1418 : {
1419 54 : eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
1420 : GF_Read, 0, iSrcLine, m_nSplitAndSwap, 1,
1421 : reinterpret_cast<void *>(reinterpret_cast<GByte *>(pData) +
1422 108 : (m_nXSize - m_nSplitAndSwap) *
1423 54 : GDALGetDataTypeSizeBytes(eReqDT)),
1424 : m_nSplitAndSwap, 1, eReqDT, 0, 0, nullptr);
1425 54 : if (eErr != CE_None)
1426 : {
1427 0 : CPLFree(pData);
1428 0 : GDALDestroyScaledProgress(pScaledProgressData);
1429 0 : return false;
1430 : }
1431 : }
1432 456 : if (m_fValOffset != 0.0)
1433 : {
1434 6 : if (eReqDT == GDT_Float32)
1435 : {
1436 12 : for (int j = 0; j < m_nXSize; j++)
1437 : {
1438 8 : static_cast<float *>(pData)[j] += m_fValOffset;
1439 : }
1440 : }
1441 : else
1442 : {
1443 6 : for (int j = 0; j < m_nXSize; j++)
1444 : {
1445 4 : static_cast<double *>(pData)[j] += m_fValOffset;
1446 : }
1447 : }
1448 : }
1449 : #ifdef CPL_LSB
1450 456 : GDALSwapWords(pData, GDALGetDataTypeSizeBytes(eReqDT), m_nXSize,
1451 : GDALGetDataTypeSizeBytes(eReqDT));
1452 : #endif
1453 456 : if (VSIFWriteL(pData, 1, nBufferSize, m_fp) != nBufferSize)
1454 : {
1455 0 : CPLFree(pData);
1456 0 : GDALDestroyScaledProgress(pScaledProgressData);
1457 0 : return false;
1458 : }
1459 456 : if (!GDALScaledProgress(static_cast<double>(i + 1) / m_nYSize, nullptr,
1460 : pScaledProgressData))
1461 : {
1462 0 : CPLFree(pData);
1463 0 : GDALDestroyScaledProgress(pScaledProgressData);
1464 0 : return false;
1465 : }
1466 : }
1467 30 : GDALDestroyScaledProgress(pScaledProgressData);
1468 30 : CPLFree(pData);
1469 :
1470 30 : return true;
1471 : }
1472 :
1473 : /************************************************************************/
1474 : /* WrapArrayAsMemDataset() */
1475 : /************************************************************************/
1476 :
1477 28 : static GDALDataset *WrapArrayAsMemDataset(int nXSize, int nYSize,
1478 : GDALDataType eReducedDT, void *pData)
1479 : {
1480 28 : CPLAssert(eReducedDT == GDT_Byte || eReducedDT == GDT_UInt16);
1481 : auto poMEMDS =
1482 28 : MEMDataset::Create("", nXSize, nYSize, 0, eReducedDT, nullptr);
1483 28 : GByte *pabyData = static_cast<GByte *>(pData);
1484 : #ifdef CPL_MSB
1485 : if (eReducedDT == GDT_Byte)
1486 : pabyData++;
1487 : #endif
1488 : auto hBand =
1489 28 : MEMCreateRasterBandEx(poMEMDS, 1, pabyData, eReducedDT, 2, 0, false);
1490 28 : poMEMDS->AddMEMBand(hBand);
1491 28 : return poMEMDS;
1492 : }
1493 :
1494 : /************************************************************************/
1495 : /* GetRoundedToUpperPowerOfTwo() */
1496 : /************************************************************************/
1497 :
1498 13 : static int GetRoundedToUpperPowerOfTwo(int nBits)
1499 : {
1500 13 : if (nBits == 3)
1501 1 : nBits = 4;
1502 12 : else if (nBits > 4 && nBits < 8)
1503 2 : nBits = 8;
1504 10 : else if (nBits > 8 && nBits < 15)
1505 1 : nBits = 16;
1506 13 : return nBits;
1507 : }
1508 :
1509 : /************************************************************************/
1510 : /* GetScaledData() */
1511 : /************************************************************************/
1512 :
1513 28 : static GUInt16 *GetScaledData(GUInt32 nDataPoints, const float *pafData,
1514 : float fMin, float fMax, double dfDecimalScale,
1515 : double dfMinScaled,
1516 : bool bOnlyPowerOfTwoDepthAllowed, int &nBits,
1517 : GInt16 &nBinaryScaleFactor)
1518 : {
1519 28 : bool bDone = false;
1520 28 : nBinaryScaleFactor = 0;
1521 : GUInt16 *panData = static_cast<GUInt16 *>(
1522 28 : VSI_MALLOC2_VERBOSE(nDataPoints, sizeof(GUInt16)));
1523 28 : if (panData == nullptr)
1524 : {
1525 0 : return nullptr;
1526 : }
1527 :
1528 28 : const double dfScaledMaxDiff = (fMax - fMin) * dfDecimalScale;
1529 28 : if (nBits == 0)
1530 : {
1531 12 : nBits = (g2int)ceil(log(ceil(dfScaledMaxDiff)) / log(2.0));
1532 12 : if (nBits > 16)
1533 : {
1534 2 : CPLError(CE_Warning, CPLE_AppDefined,
1535 : "More than 16 bits of integer precision would be "
1536 : "required. Dropping precision to fit on 16 bits");
1537 2 : nBits = 16;
1538 : }
1539 : else
1540 : {
1541 10 : bDone = true;
1542 527498 : for (GUInt32 i = 0; i < nDataPoints; i++)
1543 : {
1544 527488 : panData[i] = static_cast<GUInt16>(
1545 527488 : 0.5 + (pafData[i] * dfDecimalScale - dfMinScaled));
1546 : }
1547 : }
1548 : }
1549 :
1550 28 : if (bOnlyPowerOfTwoDepthAllowed)
1551 13 : nBits = GetRoundedToUpperPowerOfTwo(nBits);
1552 :
1553 28 : if (!bDone && nBits != 0)
1554 : {
1555 18 : if (nBits > 16)
1556 : {
1557 0 : CPLError(CE_Warning, CPLE_AppDefined,
1558 : "Maximum bit depth supported is 16. Using that");
1559 0 : nBits = 16;
1560 : }
1561 18 : const int nMaxNum = (1 << nBits) - 1;
1562 18 : double dfTemp = log(nMaxNum / dfScaledMaxDiff) / log(2.0);
1563 18 : nBinaryScaleFactor = static_cast<GInt16>(ceil(-dfTemp));
1564 18 : double dfBinaryScale = pow(2.0, -1.0 * nBinaryScaleFactor);
1565 5634 : for (GUInt32 i = 0; i < nDataPoints; i++)
1566 : {
1567 5616 : panData[i] = static_cast<GUInt16>(
1568 5616 : 0.5 +
1569 5616 : (pafData[i] * dfDecimalScale - dfMinScaled) * dfBinaryScale);
1570 : }
1571 : }
1572 :
1573 28 : return panData;
1574 : }
1575 :
1576 : /************************************************************************/
1577 : /* WritePNG() */
1578 : /************************************************************************/
1579 :
1580 : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-41.shtml
1581 14 : bool GRIB2Section567Writer::WritePNG()
1582 : {
1583 14 : float *pafData = GetFloatData();
1584 14 : if (pafData == nullptr)
1585 0 : return false;
1586 :
1587 14 : if (m_bUseZeroBits)
1588 : {
1589 : // Section 5: Data Representation Section
1590 1 : WriteUInt32(m_fp, 21); // section size
1591 1 : WriteByte(m_fp, 5); // section number
1592 1 : WriteUInt32(m_fp, m_nDataPoints);
1593 1 : WriteUInt16(m_fp, GS5_PNG);
1594 1 : WriteFloat32(
1595 : m_fp,
1596 1 : static_cast<float>(m_dfMinScaled / m_dfDecimalScale)); // ref value
1597 1 : WriteInt16(m_fp, 0); // Binary scale factor (E)
1598 1 : WriteInt16(m_fp, 0); // Decimal scale factor (D)
1599 1 : WriteByte(m_fp, 0); // Number of bits
1600 : // Type of original data: 0=Floating, 1=Integer
1601 1 : WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
1602 :
1603 : // Section 6: Bitmap section
1604 1 : WriteUInt32(m_fp, 6); // section size
1605 1 : WriteByte(m_fp, 6); // section number
1606 1 : WriteByte(m_fp, GRIB2MISSING_u1); // no bitmap
1607 :
1608 : // Section 7: Data Section
1609 1 : WriteUInt32(m_fp, 5); // section size
1610 1 : WriteByte(m_fp, 7); // section number
1611 :
1612 1 : CPLFree(pafData);
1613 :
1614 1 : return true;
1615 : }
1616 :
1617 : GDALDriver *poPNGDriver =
1618 13 : reinterpret_cast<GDALDriver *>(GDALGetDriverByName("PNG"));
1619 13 : if (poPNGDriver == nullptr)
1620 : {
1621 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find PNG driver");
1622 0 : return false;
1623 : }
1624 :
1625 13 : GInt16 nBinaryScaleFactor = 0;
1626 : GUInt16 *panData =
1627 26 : GetScaledData(m_nDataPoints, pafData, m_fMin, m_fMax, m_dfDecimalScale,
1628 13 : m_dfMinScaled, true, m_nBits, nBinaryScaleFactor);
1629 13 : if (panData == nullptr)
1630 : {
1631 0 : VSIFree(pafData);
1632 0 : return false;
1633 : }
1634 :
1635 13 : CPLFree(pafData);
1636 :
1637 26 : CPLStringList aosPNGOptions;
1638 13 : aosPNGOptions.SetNameValue("NBITS", CPLSPrintf("%d", m_nBits));
1639 :
1640 13 : const GDALDataType eReducedDT = (m_nBits <= 8) ? GDT_Byte : GDT_UInt16;
1641 : GDALDataset *poMEMDS =
1642 13 : WrapArrayAsMemDataset(m_nXSize, m_nYSize, eReducedDT, panData);
1643 :
1644 26 : const CPLString osTmpFile(VSIMemGenerateHiddenFilename("grib_driver.png"));
1645 13 : GDALDataset *poPNGDS = poPNGDriver->CreateCopy(
1646 13 : osTmpFile, poMEMDS, FALSE, aosPNGOptions.List(), nullptr, nullptr);
1647 13 : if (poPNGDS == nullptr)
1648 : {
1649 0 : CPLError(CE_Failure, CPLE_AppDefined, "PNG compression failed");
1650 0 : VSIUnlink(osTmpFile);
1651 0 : delete poMEMDS;
1652 0 : CPLFree(panData);
1653 0 : return false;
1654 : }
1655 13 : delete poPNGDS;
1656 13 : delete poMEMDS;
1657 13 : CPLFree(panData);
1658 :
1659 : // Section 5: Data Representation Section
1660 13 : WriteUInt32(m_fp, 21); // section size
1661 13 : WriteByte(m_fp, 5); // section number
1662 13 : WriteUInt32(m_fp, m_nDataPoints);
1663 13 : WriteUInt16(m_fp, GS5_PNG);
1664 13 : WriteFloat32(m_fp, static_cast<float>(m_dfMinScaled));
1665 13 : WriteInt16(m_fp, nBinaryScaleFactor); // Binary scale factor (E)
1666 13 : WriteInt16(m_fp, m_nDecimalScaleFactor); // Decimal scale factor (D)
1667 13 : WriteByte(m_fp, m_nBits); // Number of bits
1668 : // Type of original data: 0=Floating, 1=Integer
1669 13 : WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
1670 :
1671 : // Section 6: Bitmap section
1672 13 : WriteUInt32(m_fp, 6); // section size
1673 13 : WriteByte(m_fp, 6); // section number
1674 13 : WriteByte(m_fp, GRIB2MISSING_u1); // no bitmap
1675 :
1676 : // Section 7: Data Section
1677 13 : vsi_l_offset nDataLength = 0;
1678 13 : GByte *pabyData = VSIGetMemFileBuffer(osTmpFile, &nDataLength, FALSE);
1679 13 : WriteUInt32(m_fp, static_cast<GUInt32>(5 + nDataLength)); // section size
1680 13 : WriteByte(m_fp, 7); // section number
1681 13 : const size_t nDataLengthSize = static_cast<size_t>(nDataLength);
1682 : const bool bOK =
1683 13 : VSIFWriteL(pabyData, 1, nDataLengthSize, m_fp) == nDataLengthSize;
1684 :
1685 13 : VSIUnlink(osTmpFile);
1686 13 : VSIUnlink((osTmpFile + ".aux.xml").c_str());
1687 :
1688 13 : return bOK;
1689 : }
1690 :
1691 : /************************************************************************/
1692 : /* WriteJPEG2000() */
1693 : /************************************************************************/
1694 :
1695 : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp5-40.shtml
1696 17 : bool GRIB2Section567Writer::WriteJPEG2000(char **papszOptions)
1697 : {
1698 17 : float *pafData = GetFloatData();
1699 17 : if (pafData == nullptr)
1700 0 : return false;
1701 :
1702 17 : if (m_bUseZeroBits)
1703 : {
1704 : // Section 5: Data Representation Section
1705 1 : WriteUInt32(m_fp, 23); // section size
1706 1 : WriteByte(m_fp, 5); // section number
1707 1 : WriteUInt32(m_fp, m_nDataPoints);
1708 1 : WriteUInt16(m_fp, GS5_JPEG2000);
1709 1 : WriteFloat32(
1710 : m_fp,
1711 1 : static_cast<float>(m_dfMinScaled / m_dfDecimalScale)); // ref val
1712 1 : WriteInt16(m_fp, 0); // Binary scale factor (E)
1713 1 : WriteInt16(m_fp, 0); // Decimal scale factor (D)
1714 1 : WriteByte(m_fp, 0); // Number of bits
1715 : // Type of original data: 0=Floating, 1=Integer
1716 1 : WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
1717 1 : WriteByte(m_fp, 0); // compression type: lossless
1718 1 : WriteByte(m_fp, GRIB2MISSING_u1); // compression ratio
1719 :
1720 : // Section 6: Bitmap section
1721 1 : WriteUInt32(m_fp, 6); // section size
1722 1 : WriteByte(m_fp, 6); // section number
1723 1 : WriteByte(m_fp, GRIB2MISSING_u1); // no bitmap
1724 :
1725 : // Section 7: Data Section
1726 1 : WriteUInt32(m_fp, 5); // section size
1727 1 : WriteByte(m_fp, 7); // section number
1728 :
1729 1 : CPLFree(pafData);
1730 :
1731 1 : return true;
1732 : }
1733 :
1734 16 : GDALDriver *poJ2KDriver = nullptr;
1735 16 : const char *pszJ2KDriver = GetBandOption(papszOptions, nullptr, m_nBand,
1736 : "JPEG2000_DRIVER", nullptr);
1737 16 : if (pszJ2KDriver)
1738 : {
1739 : poJ2KDriver =
1740 6 : reinterpret_cast<GDALDriver *>(GDALGetDriverByName(pszJ2KDriver));
1741 : }
1742 : else
1743 : {
1744 20 : for (size_t i = 0; i < CPL_ARRAYSIZE(apszJ2KDrivers); i++)
1745 : {
1746 : poJ2KDriver = reinterpret_cast<GDALDriver *>(
1747 20 : GDALGetDriverByName(apszJ2KDrivers[i]));
1748 20 : if (poJ2KDriver)
1749 : {
1750 10 : CPLDebug("GRIB", "Using %s", poJ2KDriver->GetDescription());
1751 10 : break;
1752 : }
1753 : }
1754 : }
1755 16 : if (poJ2KDriver == nullptr)
1756 : {
1757 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find JPEG2000 driver");
1758 1 : VSIFree(pafData);
1759 1 : return false;
1760 : }
1761 :
1762 15 : GInt16 nBinaryScaleFactor = 0;
1763 : GUInt16 *panData =
1764 30 : GetScaledData(m_nDataPoints, pafData, m_fMin, m_fMax, m_dfDecimalScale,
1765 15 : m_dfMinScaled, false, m_nBits, nBinaryScaleFactor);
1766 15 : if (panData == nullptr)
1767 : {
1768 0 : VSIFree(pafData);
1769 0 : return false;
1770 : }
1771 :
1772 15 : CPLFree(pafData);
1773 :
1774 30 : CPLStringList aosJ2KOptions;
1775 15 : int nCompressionRatio = atoi(GetBandOption(papszOptions, nullptr, m_nBand,
1776 : "COMPRESSION_RATIO", "1"));
1777 15 : if (m_nDataPoints < 10000 && nCompressionRatio > 1)
1778 : {
1779 : // Lossy compression with too few pixels is really lossy due to how
1780 : // codec work
1781 1 : CPLDebug("GRIB", "Forcing JPEG2000 lossless mode given "
1782 : "the low number of pixels");
1783 1 : nCompressionRatio = 1;
1784 : }
1785 15 : const bool bLossLess = nCompressionRatio <= 1;
1786 15 : if (EQUAL(poJ2KDriver->GetDescription(), "JP2KAK"))
1787 : {
1788 0 : if (bLossLess)
1789 : {
1790 0 : aosJ2KOptions.SetNameValue("QUALITY", "100");
1791 : }
1792 : else
1793 : {
1794 : aosJ2KOptions.SetNameValue(
1795 : "QUALITY",
1796 0 : CPLSPrintf("%d", std::max(1, 100 / nCompressionRatio)));
1797 : }
1798 : }
1799 15 : else if (EQUAL(poJ2KDriver->GetDescription(), "JP2OPENJPEG"))
1800 : {
1801 12 : if (bLossLess)
1802 : {
1803 11 : aosJ2KOptions.SetNameValue("QUALITY", "100");
1804 11 : aosJ2KOptions.SetNameValue("REVERSIBLE", "YES");
1805 : }
1806 : else
1807 : {
1808 : aosJ2KOptions.SetNameValue(
1809 1 : "QUALITY", CPLSPrintf("%f", 100.0 / nCompressionRatio));
1810 : }
1811 : }
1812 3 : else if (EQUAL(poJ2KDriver->GetDescription(), "JP2ECW"))
1813 : {
1814 2 : if (bLossLess)
1815 : {
1816 1 : aosJ2KOptions.SetNameValue("TARGET", "0");
1817 : }
1818 : else
1819 : {
1820 : aosJ2KOptions.SetNameValue(
1821 1 : "TARGET", CPLSPrintf("%f", 100.0 - 100.0 / nCompressionRatio));
1822 : }
1823 : }
1824 15 : aosJ2KOptions.SetNameValue("NBITS", CPLSPrintf("%d", m_nBits));
1825 :
1826 15 : const GDALDataType eReducedDT = (m_nBits <= 8) ? GDT_Byte : GDT_UInt16;
1827 : GDALDataset *poMEMDS =
1828 15 : WrapArrayAsMemDataset(m_nXSize, m_nYSize, eReducedDT, panData);
1829 :
1830 30 : const CPLString osTmpFile(VSIMemGenerateHiddenFilename("grib_driver.j2k"));
1831 15 : GDALDataset *poJ2KDS = poJ2KDriver->CreateCopy(
1832 15 : osTmpFile, poMEMDS, FALSE, aosJ2KOptions.List(), nullptr, nullptr);
1833 15 : if (poJ2KDS == nullptr)
1834 : {
1835 1 : CPLError(CE_Failure, CPLE_AppDefined, "JPEG2000 compression failed");
1836 1 : VSIUnlink(osTmpFile);
1837 1 : delete poMEMDS;
1838 1 : CPLFree(panData);
1839 1 : return false;
1840 : }
1841 14 : delete poJ2KDS;
1842 14 : delete poMEMDS;
1843 14 : CPLFree(panData);
1844 :
1845 : // Section 5: Data Representation Section
1846 14 : WriteUInt32(m_fp, 23); // section size
1847 14 : WriteByte(m_fp, 5); // section number
1848 14 : WriteUInt32(m_fp, m_nDataPoints);
1849 14 : WriteUInt16(m_fp, GS5_JPEG2000);
1850 14 : WriteFloat32(m_fp, static_cast<float>(m_dfMinScaled));
1851 14 : WriteInt16(m_fp, nBinaryScaleFactor); // Binary scale factor (E)
1852 14 : WriteInt16(m_fp, m_nDecimalScaleFactor); // Decimal scale factor (D)
1853 14 : WriteByte(m_fp, m_nBits); // Number of bits
1854 : // Type of original data: 0=Floating, 1=Integer
1855 14 : WriteByte(m_fp, GDALDataTypeIsFloating(m_eDT) ? 0 : 1);
1856 : // compression type: lossless(0) or lossy(1)
1857 14 : WriteByte(m_fp, bLossLess ? 0 : 1);
1858 14 : WriteByte(m_fp, bLossLess ? GRIB2MISSING_u1
1859 : : nCompressionRatio); // compression ratio
1860 :
1861 : // Section 6: Bitmap section
1862 14 : WriteUInt32(m_fp, 6); // section size
1863 14 : WriteByte(m_fp, 6); // section number
1864 14 : WriteByte(m_fp, GRIB2MISSING_u1); // no bitmap
1865 :
1866 : // Section 7: Data Section
1867 14 : vsi_l_offset nDataLength = 0;
1868 14 : GByte *pabyData = VSIGetMemFileBuffer(osTmpFile, &nDataLength, FALSE);
1869 14 : WriteUInt32(m_fp, static_cast<GUInt32>(5 + nDataLength)); // section size
1870 14 : WriteByte(m_fp, 7); // section number
1871 14 : const size_t nDataLengthSize = static_cast<size_t>(nDataLength);
1872 : const bool bOK =
1873 14 : VSIFWriteL(pabyData, 1, nDataLengthSize, m_fp) == nDataLengthSize;
1874 :
1875 14 : VSIUnlink(osTmpFile);
1876 14 : VSIUnlink((osTmpFile + ".aux.xml").c_str());
1877 :
1878 14 : return bOK;
1879 : }
1880 :
1881 : /************************************************************************/
1882 : /* Write() */
1883 : /************************************************************************/
1884 :
1885 156 : bool GRIB2Section567Writer::Write(float fValOffset, char **papszOptions,
1886 : GDALProgressFunc pfnProgress,
1887 : void *pProgressData)
1888 : {
1889 156 : m_fValOffset = fValOffset;
1890 :
1891 : typedef enum
1892 : {
1893 : SIMPLE_PACKING,
1894 : COMPLEX_PACKING,
1895 : IEEE_FLOATING_POINT,
1896 : PNG,
1897 : JPEG2000
1898 : } GRIBDataEncoding;
1899 :
1900 156 : if (m_eDT != GDT_Byte && m_eDT != GDT_UInt16 && m_eDT != GDT_Int16 &&
1901 56 : m_eDT != GDT_UInt32 && m_eDT != GDT_Int32 && m_eDT != GDT_Float32 &&
1902 30 : m_eDT != GDT_Float64)
1903 : {
1904 5 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported data type: %s",
1905 : GDALGetDataTypeName(m_eDT));
1906 5 : return false;
1907 : }
1908 : const char *pszDataEncoding =
1909 151 : GetBandOption(papszOptions, nullptr, m_nBand, "DATA_ENCODING", "AUTO");
1910 151 : GRIBDataEncoding eDataEncoding(SIMPLE_PACKING);
1911 151 : const char *pszJ2KDriver = GetBandOption(papszOptions, nullptr, m_nBand,
1912 : "JPEG2000_DRIVER", nullptr);
1913 151 : const char *pszSpatialDifferencingOrder = GetBandOption(
1914 : papszOptions, nullptr, m_nBand, "SPATIAL_DIFFERENCING_ORDER", nullptr);
1915 151 : if (pszJ2KDriver && pszSpatialDifferencingOrder)
1916 : {
1917 1 : CPLError(CE_Failure, CPLE_NotSupported,
1918 : "JPEG2000_DRIVER and SPATIAL_DIFFERENCING_ORDER are not "
1919 : "compatible");
1920 1 : return false;
1921 : }
1922 :
1923 150 : if (m_bHasNoData && !EQUAL(pszDataEncoding, "COMPLEX_PACKING") &&
1924 : pszSpatialDifferencingOrder == nullptr)
1925 : {
1926 : double *padfVals = static_cast<double *>(
1927 6 : VSI_MALLOC2_VERBOSE(m_nXSize, sizeof(double)));
1928 6 : if (padfVals == nullptr)
1929 0 : return false;
1930 6 : bool bFoundNoData = false;
1931 7 : for (int j = 0; j < m_nYSize; j++)
1932 : {
1933 6 : CPLErr eErr = m_poSrcDS->GetRasterBand(m_nBand)->RasterIO(
1934 : GF_Read, 0, j, m_nXSize, 1, padfVals, m_nXSize, 1, GDT_Float64,
1935 : 0, 0, nullptr);
1936 6 : if (eErr != CE_None)
1937 : {
1938 0 : VSIFree(padfVals);
1939 0 : return false;
1940 : }
1941 7 : for (int i = 0; i < m_nXSize; i++)
1942 : {
1943 6 : if (padfVals[i] == m_dfNoData)
1944 : {
1945 5 : bFoundNoData = true;
1946 5 : break;
1947 : }
1948 : }
1949 6 : if (bFoundNoData)
1950 5 : break;
1951 : }
1952 6 : VSIFree(padfVals);
1953 :
1954 6 : if (!bFoundNoData)
1955 : {
1956 1 : m_bHasNoData = false;
1957 : }
1958 : }
1959 :
1960 150 : if (EQUAL(pszDataEncoding, "AUTO"))
1961 : {
1962 81 : if (m_bHasNoData || pszSpatialDifferencingOrder != nullptr)
1963 : {
1964 5 : eDataEncoding = COMPLEX_PACKING;
1965 5 : CPLDebug("GRIB", "Using COMPLEX_PACKING");
1966 : }
1967 76 : else if (pszJ2KDriver != nullptr)
1968 : {
1969 6 : eDataEncoding = JPEG2000;
1970 6 : CPLDebug("GRIB", "Using JPEG2000");
1971 : }
1972 70 : else if (m_eDT == GDT_Float32 || m_eDT == GDT_Float64)
1973 : {
1974 20 : eDataEncoding = IEEE_FLOATING_POINT;
1975 20 : CPLDebug("GRIB", "Using IEEE_FLOATING_POINT");
1976 : }
1977 : else
1978 : {
1979 50 : CPLDebug("GRIB", "Using SIMPLE_PACKING");
1980 : }
1981 : }
1982 69 : else if (EQUAL(pszDataEncoding, "SIMPLE_PACKING"))
1983 : {
1984 16 : eDataEncoding = SIMPLE_PACKING;
1985 : }
1986 53 : else if (EQUAL(pszDataEncoding, "COMPLEX_PACKING"))
1987 : {
1988 17 : eDataEncoding = COMPLEX_PACKING;
1989 : }
1990 36 : else if (EQUAL(pszDataEncoding, "IEEE_FLOATING_POINT"))
1991 : {
1992 10 : eDataEncoding = IEEE_FLOATING_POINT;
1993 : }
1994 26 : else if (EQUAL(pszDataEncoding, "PNG"))
1995 : {
1996 14 : eDataEncoding = PNG;
1997 : }
1998 12 : else if (EQUAL(pszDataEncoding, "JPEG2000"))
1999 : {
2000 11 : eDataEncoding = JPEG2000;
2001 : }
2002 : else
2003 : {
2004 1 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported DATA_ENCODING=%s",
2005 : pszDataEncoding);
2006 1 : return false;
2007 : }
2008 :
2009 : const char *pszBits =
2010 149 : GetBandOption(papszOptions, nullptr, m_nBand, "NBITS", nullptr);
2011 149 : if (pszBits == nullptr && eDataEncoding != IEEE_FLOATING_POINT)
2012 : {
2013 200 : pszBits = m_poSrcDS->GetRasterBand(m_nBand)->GetMetadataItem(
2014 100 : "DRS_NBITS", "GRIB");
2015 : }
2016 49 : else if (pszBits != nullptr && eDataEncoding == IEEE_FLOATING_POINT)
2017 : {
2018 1 : CPLError(CE_Warning, CPLE_NotSupported,
2019 : "NBITS ignored for DATA_ENCODING = IEEE_FLOATING_POINT");
2020 : }
2021 149 : if (pszBits == nullptr)
2022 : {
2023 123 : pszBits = "0";
2024 : }
2025 149 : m_nBits = std::max(0, atoi(pszBits));
2026 149 : if (m_nBits > 31)
2027 : {
2028 1 : CPLError(CE_Warning, CPLE_NotSupported, "NBITS clamped to 31");
2029 1 : m_nBits = 31;
2030 : }
2031 :
2032 149 : const char *pszDecimalScaleFactor = GetBandOption(
2033 : papszOptions, nullptr, m_nBand, "DECIMAL_SCALE_FACTOR", nullptr);
2034 149 : if (pszDecimalScaleFactor != nullptr)
2035 : {
2036 10 : m_nDecimalScaleFactor = atoi(pszDecimalScaleFactor);
2037 10 : if (m_nDecimalScaleFactor != 0 && eDataEncoding == IEEE_FLOATING_POINT)
2038 : {
2039 1 : CPLError(CE_Warning, CPLE_NotSupported,
2040 : "DECIMAL_SCALE_FACTOR ignored for "
2041 : "DATA_ENCODING = IEEE_FLOATING_POINT");
2042 : }
2043 9 : else if (m_nDecimalScaleFactor > 0 && !GDALDataTypeIsFloating(m_eDT))
2044 : {
2045 1 : CPLError(CE_Warning, CPLE_AppDefined,
2046 : "DECIMAL_SCALE_FACTOR > 0 makes no sense for integer "
2047 : "data types. Ignored");
2048 1 : m_nDecimalScaleFactor = 0;
2049 : }
2050 : }
2051 139 : else if (eDataEncoding != IEEE_FLOATING_POINT)
2052 : {
2053 : pszDecimalScaleFactor =
2054 110 : m_poSrcDS->GetRasterBand(m_nBand)->GetMetadataItem(
2055 110 : "DRS_DECIMAL_SCALE_FACTOR", "GRIB");
2056 110 : if (pszDecimalScaleFactor != nullptr)
2057 : {
2058 6 : m_nDecimalScaleFactor = atoi(pszDecimalScaleFactor);
2059 : }
2060 : }
2061 149 : m_dfDecimalScale = pow(10.0, static_cast<double>(m_nDecimalScaleFactor));
2062 :
2063 149 : if (pszJ2KDriver != nullptr && eDataEncoding != JPEG2000)
2064 : {
2065 2 : CPLError(CE_Warning, CPLE_AppDefined,
2066 : "JPEG2000_DRIVER option ignored for "
2067 : "non-JPEG2000 DATA_ENCODING");
2068 : }
2069 149 : if (pszSpatialDifferencingOrder && eDataEncoding != COMPLEX_PACKING)
2070 : {
2071 1 : CPLError(CE_Warning, CPLE_AppDefined,
2072 : "SPATIAL_DIFFERENCING_ORDER option ignored for "
2073 : "non-COMPLEX_PACKING DATA_ENCODING");
2074 : }
2075 149 : if (m_bHasNoData && eDataEncoding != COMPLEX_PACKING)
2076 : {
2077 2 : CPLError(CE_Warning, CPLE_AppDefined,
2078 : "non-COMPLEX_PACKING DATA_ENCODING cannot preserve nodata");
2079 : }
2080 :
2081 149 : if (eDataEncoding == SIMPLE_PACKING)
2082 : {
2083 66 : return WriteSimplePacking();
2084 : }
2085 83 : else if (eDataEncoding == COMPLEX_PACKING)
2086 : {
2087 22 : const int nSpatialDifferencingOrder =
2088 22 : pszSpatialDifferencingOrder ? atoi(pszSpatialDifferencingOrder) : 0;
2089 22 : return WriteComplexPacking(nSpatialDifferencingOrder);
2090 : }
2091 61 : else if (eDataEncoding == IEEE_FLOATING_POINT)
2092 : {
2093 30 : return WriteIEEE(pfnProgress, pProgressData);
2094 : }
2095 31 : else if (eDataEncoding == PNG)
2096 : {
2097 14 : return WritePNG();
2098 : }
2099 : else /* if( eDataEncoding == JPEG2000 ) */
2100 : {
2101 17 : return WriteJPEG2000(papszOptions);
2102 : }
2103 : }
2104 :
2105 : /************************************************************************/
2106 : /* GetIDSOption() */
2107 : /************************************************************************/
2108 :
2109 1134 : static const char *GetIDSOption(char **papszOptions, GDALDataset *poSrcDS,
2110 : int nBand, const char *pszKey,
2111 : const char *pszDefault)
2112 : {
2113 : const char *pszValue =
2114 1134 : GetBandOption(papszOptions, nullptr, nBand,
2115 2268 : (CPLString("IDS_") + pszKey).c_str(), nullptr);
2116 1134 : if (pszValue == nullptr)
2117 : {
2118 : const char *pszIDS =
2119 1133 : GetBandOption(papszOptions, poSrcDS, nBand, "IDS", nullptr);
2120 1133 : if (pszIDS != nullptr)
2121 : {
2122 153 : char **papszTokens = CSLTokenizeString2(pszIDS, " ", 0);
2123 153 : pszValue = CSLFetchNameValue(papszTokens, pszKey);
2124 153 : if (pszValue)
2125 142 : pszValue = CPLSPrintf("%s", pszValue);
2126 153 : CSLDestroy(papszTokens);
2127 : }
2128 : }
2129 1134 : if (pszValue == nullptr)
2130 991 : pszValue = pszDefault;
2131 1134 : return pszValue;
2132 : }
2133 :
2134 : /************************************************************************/
2135 : /* WriteSection1() */
2136 : /************************************************************************/
2137 :
2138 162 : static void WriteSection1(VSILFILE *fp, GDALDataset *poSrcDS, int nBand,
2139 : char **papszOptions)
2140 : {
2141 : // Section 1: Identification Section
2142 162 : WriteUInt32(fp, 21); // section size
2143 162 : WriteByte(fp, 1); // section number
2144 :
2145 : GUInt16 nCenter = static_cast<GUInt16>(
2146 162 : atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "CENTER",
2147 162 : CPLSPrintf("%d", GRIB2MISSING_u1))));
2148 162 : WriteUInt16(fp, nCenter);
2149 :
2150 : GUInt16 nSubCenter = static_cast<GUInt16>(
2151 162 : atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "SUBCENTER",
2152 162 : CPLSPrintf("%d", GRIB2MISSING_u2))));
2153 162 : WriteUInt16(fp, nSubCenter);
2154 :
2155 : GByte nMasterTable = static_cast<GByte>(
2156 162 : atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "MASTER_TABLE", "2")));
2157 162 : WriteByte(fp, nMasterTable);
2158 :
2159 162 : WriteByte(fp, 0); // local table
2160 :
2161 162 : GByte nSignfRefTime = static_cast<GByte>(atoi(
2162 162 : GetIDSOption(papszOptions, poSrcDS, nBand, "SIGNF_REF_TIME", "0")));
2163 162 : WriteByte(fp, nSignfRefTime); // Significance of reference time
2164 :
2165 : const char *pszRefTime =
2166 162 : GetIDSOption(papszOptions, poSrcDS, nBand, "REF_TIME", "");
2167 162 : int nYear = 1970, nMonth = 1, nDay = 1, nHour = 0, nMinute = 0, nSecond = 0;
2168 162 : sscanf(pszRefTime, "%04d-%02d-%02dT%02d:%02d:%02dZ", &nYear, &nMonth, &nDay,
2169 : &nHour, &nMinute, &nSecond);
2170 162 : WriteUInt16(fp, nYear);
2171 162 : WriteByte(fp, nMonth);
2172 162 : WriteByte(fp, nDay);
2173 162 : WriteByte(fp, nHour);
2174 162 : WriteByte(fp, nMinute);
2175 162 : WriteByte(fp, nSecond);
2176 :
2177 : GByte nProdStatus = static_cast<GByte>(
2178 162 : atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "PROD_STATUS",
2179 162 : CPLSPrintf("%d", GRIB2MISSING_u1))));
2180 162 : WriteByte(fp, nProdStatus);
2181 :
2182 : GByte nType = static_cast<GByte>(
2183 162 : atoi(GetIDSOption(papszOptions, poSrcDS, nBand, "TYPE",
2184 162 : CPLSPrintf("%d", GRIB2MISSING_u1))));
2185 162 : WriteByte(fp, nType);
2186 162 : }
2187 :
2188 : /************************************************************************/
2189 : /* WriteAssembledPDS() */
2190 : /************************************************************************/
2191 :
2192 12 : static void WriteAssembledPDS(VSILFILE *fp, const gtemplate *mappds,
2193 : bool bWriteExt, char **papszTokens,
2194 : std::vector<int> &anVals)
2195 : {
2196 12 : const int iStart = bWriteExt ? mappds->maplen : 0;
2197 12 : const int iEnd =
2198 12 : bWriteExt ? mappds->maplen + mappds->extlen : mappds->maplen;
2199 173 : for (int i = iStart; i < iEnd; i++)
2200 : {
2201 161 : const int nVal = atoi(papszTokens[i]);
2202 161 : anVals.push_back(nVal);
2203 161 : const int nEltSize =
2204 161 : bWriteExt ? mappds->ext[i - mappds->maplen] : mappds->map[i];
2205 161 : if (nEltSize == 1)
2206 : {
2207 90 : if (nVal < 0 || nVal > 255)
2208 : {
2209 2 : CPLError(CE_Warning, CPLE_AppDefined,
2210 : "Value %d of index %d in PDS should be in [0,255] "
2211 : "range",
2212 : nVal, i);
2213 : }
2214 90 : WriteByte(fp, nVal);
2215 : }
2216 71 : else if (nEltSize == 2)
2217 : {
2218 25 : if (nVal < 0 || nVal > 65535)
2219 : {
2220 2 : CPLError(CE_Warning, CPLE_AppDefined,
2221 : "Value %d of index %d in PDS should be in [0,65535] "
2222 : "range",
2223 : nVal, i);
2224 : }
2225 25 : WriteUInt16(fp, nVal);
2226 : }
2227 46 : else if (nEltSize == 4)
2228 : {
2229 6 : GIntBig nBigVal = CPLAtoGIntBig(papszTokens[i]);
2230 6 : anVals.back() = static_cast<int>(nBigVal);
2231 6 : if (nBigVal < 0 || nBigVal > static_cast<GIntBig>(UINT_MAX))
2232 : {
2233 0 : CPLError(CE_Warning, CPLE_AppDefined,
2234 : "Value " CPL_FRMT_GIB " of index %d in PDS should be "
2235 : "in [0,%d] range",
2236 : nBigVal, i, INT_MAX);
2237 : }
2238 6 : WriteUInt32(fp, static_cast<GUInt32>(nBigVal));
2239 : }
2240 40 : else if (nEltSize == -1)
2241 : {
2242 16 : if (nVal < -128 || nVal > 127)
2243 : {
2244 2 : CPLError(CE_Warning, CPLE_AppDefined,
2245 : "Value %d of index %d in PDS should be in [-128,127] "
2246 : "range",
2247 : nVal, i);
2248 : }
2249 16 : WriteSByte(fp, nVal);
2250 : }
2251 24 : else if (nEltSize == -2)
2252 : {
2253 1 : if (nVal < -32768 || nVal > 32767)
2254 : {
2255 1 : CPLError(CE_Warning, CPLE_AppDefined,
2256 : "Value %d of index %d in PDS should be in "
2257 : "[-32768,32767] range",
2258 : nVal, i);
2259 : }
2260 1 : WriteInt16(fp, nVal);
2261 : }
2262 23 : else if (nEltSize == -4)
2263 : {
2264 23 : GIntBig nBigVal = CPLAtoGIntBig(papszTokens[i]);
2265 23 : if (nBigVal < INT_MIN || nBigVal > INT_MAX)
2266 : {
2267 1 : CPLError(CE_Warning, CPLE_AppDefined,
2268 : "Value " CPL_FRMT_GIB " of index %d in PDS should be "
2269 : "in [%d,%d] range",
2270 : nBigVal, i, INT_MIN, INT_MAX);
2271 : }
2272 23 : WriteInt32(fp, atoi(papszTokens[i]));
2273 : }
2274 : else
2275 : {
2276 0 : CPLAssert(false);
2277 : }
2278 : }
2279 12 : }
2280 :
2281 : /************************************************************************/
2282 : /* ComputeValOffset() */
2283 : /************************************************************************/
2284 :
2285 38 : static float ComputeValOffset(int nTokens, char **papszTokens,
2286 : const char *pszInputUnit)
2287 : {
2288 38 : float fValOffset = 0.0f;
2289 :
2290 : // Parameter category 0 = Temperature
2291 38 : if (nTokens >= 2 && atoi(papszTokens[0]) == 0)
2292 : {
2293 : // Cf
2294 : // https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-2-0-0.shtml
2295 : // PARAMETERS FOR DISCIPLINE 0 CATEGORY 0
2296 11 : int nParamNumber = atoi(papszTokens[1]);
2297 11 : if ((nParamNumber >= 0 && nParamNumber <= 18 && nParamNumber != 8 &&
2298 11 : nParamNumber != 10 && nParamNumber != 11 && nParamNumber != 16) ||
2299 1 : nParamNumber == 21 || nParamNumber == 27)
2300 : {
2301 10 : if (pszInputUnit == nullptr || EQUAL(pszInputUnit, "C") ||
2302 6 : EQUAL(pszInputUnit, "[C]"))
2303 : {
2304 8 : fValOffset = 273.15f;
2305 8 : CPLDebug("GRIB",
2306 : "Applying a %f offset to convert from "
2307 : "Celsius to Kelvin",
2308 : fValOffset);
2309 : }
2310 : }
2311 : }
2312 :
2313 38 : return fValOffset;
2314 : }
2315 :
2316 : /************************************************************************/
2317 : /* WriteSection4() */
2318 : /************************************************************************/
2319 :
2320 162 : static bool WriteSection4(VSILFILE *fp, GDALDataset *poSrcDS, int nBand,
2321 : char **papszOptions, float &fValOffset)
2322 : {
2323 : // Section 4: Product Definition Section
2324 162 : vsi_l_offset nStartSection4 = VSIFTellL(fp);
2325 162 : WriteUInt32(fp, GRIB2MISSING_u4); // section size
2326 162 : WriteByte(fp, 4); // section number
2327 162 : WriteUInt16(fp, 0); // Number of coordinate values after template
2328 :
2329 : // 0 = Analysis or forecast at a horizontal level or in a horizontal
2330 : // layer at a point in time
2331 : int nPDTN =
2332 162 : atoi(GetBandOption(papszOptions, poSrcDS, nBand, "PDS_PDTN", "0"));
2333 162 : const char *pszPDSTemplateNumbers = GetBandOption(
2334 : papszOptions, nullptr, nBand, "PDS_TEMPLATE_NUMBERS", nullptr);
2335 162 : const char *pszPDSTemplateAssembledValues = GetBandOption(
2336 : papszOptions, nullptr, nBand, "PDS_TEMPLATE_ASSEMBLED_VALUES", nullptr);
2337 162 : if (pszPDSTemplateNumbers == nullptr &&
2338 : pszPDSTemplateAssembledValues == nullptr)
2339 : {
2340 140 : pszPDSTemplateNumbers = GetBandOption(papszOptions, poSrcDS, nBand,
2341 : "PDS_TEMPLATE_NUMBERS", nullptr);
2342 : }
2343 324 : CPLString osInputUnit;
2344 : const char *pszInputUnit =
2345 162 : GetBandOption(papszOptions, nullptr, nBand, "INPUT_UNIT", nullptr);
2346 162 : if (pszInputUnit == nullptr)
2347 : {
2348 : const char *pszGribUnit =
2349 160 : poSrcDS->GetRasterBand(nBand)->GetMetadataItem("GRIB_UNIT");
2350 160 : if (pszGribUnit != nullptr)
2351 : {
2352 21 : osInputUnit = pszGribUnit;
2353 21 : pszInputUnit = osInputUnit.c_str();
2354 : }
2355 : }
2356 162 : WriteUInt16(fp, nPDTN); // PDTN
2357 162 : if (nPDTN == 0 && pszPDSTemplateNumbers == nullptr &&
2358 : pszPDSTemplateAssembledValues == nullptr)
2359 : {
2360 : // See http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_temp4-0.shtml
2361 120 : WriteByte(fp, GRIB2MISSING_u1); // Parameter category = Missing
2362 120 : WriteByte(fp, GRIB2MISSING_u1); // Parameter number = Missing
2363 120 : WriteByte(fp, GRIB2MISSING_u1); // Type of generating process = Missing
2364 120 : WriteByte(fp, 0); // Background generating process identifier
2365 : // Analysis or forecast generating process identified
2366 120 : WriteByte(fp, GRIB2MISSING_u1);
2367 120 : WriteUInt16(fp, 0); // Hours
2368 120 : WriteByte(fp, 0); // Minutes
2369 120 : WriteByte(fp, 0); // Indicator of unit of time range: 0=Minute
2370 120 : WriteUInt32(fp, 0); // Forecast time in units
2371 120 : WriteByte(fp, 0); // Type of first fixed surface
2372 120 : WriteByte(fp, 0); // Scale factor of first fixed surface
2373 120 : WriteUInt32(fp, 0); // Type of second fixed surface
2374 120 : WriteByte(fp, GRIB2MISSING_u1); // Type of second fixed surface
2375 120 : WriteByte(fp, GRIB2MISSING_u1); // Scale factor of second fixed surface
2376 : // Scaled value of second fixed surface
2377 120 : WriteUInt32(fp, GRIB2MISSING_u4);
2378 : }
2379 42 : else if (pszPDSTemplateNumbers == nullptr &&
2380 : pszPDSTemplateAssembledValues == nullptr)
2381 : {
2382 1 : CPLError(CE_Failure, CPLE_NotSupported,
2383 : "PDS_PDTN != 0 specified but both PDS_TEMPLATE_NUMBERS and "
2384 : "PDS_TEMPLATE_ASSEMBLED_VALUES missing");
2385 1 : return false;
2386 : }
2387 41 : else if (pszPDSTemplateNumbers != nullptr &&
2388 : pszPDSTemplateAssembledValues != nullptr)
2389 : {
2390 1 : CPLError(CE_Failure, CPLE_NotSupported,
2391 : "PDS_TEMPLATE_NUMBERS and "
2392 : "PDS_TEMPLATE_ASSEMBLED_VALUES are exclusive");
2393 1 : return false;
2394 : }
2395 40 : else if (pszPDSTemplateNumbers != nullptr)
2396 : {
2397 30 : char **papszTokens = CSLTokenizeString2(pszPDSTemplateNumbers, " ", 0);
2398 30 : const int nTokens = CSLCount(papszTokens);
2399 :
2400 30 : fValOffset = ComputeValOffset(nTokens, papszTokens, pszInputUnit);
2401 :
2402 952 : for (int i = 0; papszTokens[i] != nullptr; i++)
2403 : {
2404 922 : int nVal = atoi(papszTokens[i]);
2405 922 : if (nVal < 0 || nVal > 255)
2406 : {
2407 2 : CPLError(CE_Warning, CPLE_AppDefined,
2408 : "Value %d of index %d in PDS should be in [0,255] "
2409 : "range",
2410 : nVal, i);
2411 : }
2412 922 : WriteByte(fp, nVal);
2413 : }
2414 30 : CSLDestroy(papszTokens);
2415 :
2416 : // Read back section
2417 30 : PatchSectionSize(fp, nStartSection4);
2418 :
2419 30 : vsi_l_offset nCurOffset = VSIFTellL(fp);
2420 30 : VSIFSeekL(fp, nStartSection4, SEEK_SET);
2421 30 : size_t nSizeSect4 = static_cast<size_t>(nCurOffset - nStartSection4);
2422 30 : GByte *pabySect4 = static_cast<GByte *>(CPLMalloc(nSizeSect4));
2423 30 : VSIFReadL(pabySect4, 1, nSizeSect4, fp);
2424 30 : VSIFSeekL(fp, nCurOffset, SEEK_SET);
2425 :
2426 : // Check consistency with template definition
2427 30 : g2int iofst = 0;
2428 30 : g2int pdsnum = 0;
2429 30 : g2int *pdstempl = nullptr;
2430 30 : g2int mappdslen = 0;
2431 30 : g2float *coordlist = nullptr;
2432 30 : g2int numcoord = 0;
2433 : int ret =
2434 30 : g2_unpack4(pabySect4, static_cast<g2int>(nSizeSect4), &iofst,
2435 : &pdsnum, &pdstempl, &mappdslen, &coordlist, &numcoord);
2436 30 : CPLFree(pabySect4);
2437 30 : if (ret == 0)
2438 : {
2439 29 : gtemplate *mappds = extpdstemplate(pdsnum, pdstempl);
2440 29 : free(pdstempl);
2441 29 : free(coordlist);
2442 29 : if (mappds)
2443 : {
2444 29 : int nTemplateByteCount = 0;
2445 564 : for (int i = 0; i < mappds->maplen; i++)
2446 535 : nTemplateByteCount += abs(mappds->map[i]);
2447 39 : for (int i = 0; i < mappds->extlen; i++)
2448 10 : nTemplateByteCount += abs(mappds->ext[i]);
2449 29 : if (nTokens < nTemplateByteCount)
2450 : {
2451 1 : CPLError(CE_Failure, CPLE_AppDefined,
2452 : "PDS_PDTN = %d (with provided elements) requires "
2453 : "%d bytes in PDS_TEMPLATE_NUMBERS. "
2454 : "Only %d provided",
2455 : nPDTN, nTemplateByteCount, nTokens);
2456 1 : free(mappds->ext);
2457 1 : free(mappds);
2458 1 : return false;
2459 : }
2460 28 : else if (nTokens > nTemplateByteCount)
2461 : {
2462 1 : CPLError(CE_Warning, CPLE_AppDefined,
2463 : "PDS_PDTN = %d (with provided elements) requires "
2464 : "%d bytes in PDS_TEMPLATE_NUMBERS. "
2465 : "But %d provided. Extra bytes will be ignored "
2466 : "by readers",
2467 : nPDTN, nTemplateByteCount, nTokens);
2468 : }
2469 :
2470 28 : free(mappds->ext);
2471 28 : free(mappds);
2472 : }
2473 : }
2474 : else
2475 : {
2476 1 : free(pdstempl);
2477 1 : free(coordlist);
2478 1 : CPLError(CE_Warning, CPLE_AppDefined,
2479 : "PDS_PDTN = %d is unknown. Product will not be "
2480 : "correctly read by this driver (but potentially valid "
2481 : "for other readers)",
2482 : nPDTN);
2483 : }
2484 : }
2485 : else
2486 : {
2487 10 : gtemplate *mappds = getpdstemplate(nPDTN);
2488 10 : if (mappds == nullptr)
2489 : {
2490 1 : CPLError(CE_Failure, CPLE_NotSupported,
2491 : "PDS_PDTN = %d is unknown, so it is not possible to use "
2492 : "PDS_TEMPLATE_ASSEMBLED_VALUES. Use PDS_TEMPLATE_NUMBERS "
2493 : "instead",
2494 : nPDTN);
2495 3 : return false;
2496 : }
2497 : char **papszTokens =
2498 9 : CSLTokenizeString2(pszPDSTemplateAssembledValues, " ", 0);
2499 9 : const int nTokens = CSLCount(papszTokens);
2500 9 : if (nTokens < mappds->maplen)
2501 : {
2502 1 : CPLError(CE_Failure, CPLE_AppDefined,
2503 : "PDS_PDTN = %d requires at least %d elements in "
2504 : "PDS_TEMPLATE_ASSEMBLED_VALUES. Only %d provided",
2505 : nPDTN, mappds->maplen, nTokens);
2506 1 : free(mappds);
2507 1 : CSLDestroy(papszTokens);
2508 1 : return false;
2509 : }
2510 :
2511 8 : fValOffset = ComputeValOffset(nTokens, papszTokens, pszInputUnit);
2512 :
2513 8 : std::vector<int> anVals;
2514 8 : WriteAssembledPDS(fp, mappds, false, papszTokens, anVals);
2515 :
2516 8 : if (mappds->needext && !anVals.empty())
2517 : {
2518 5 : free(mappds);
2519 5 : mappds = extpdstemplate(nPDTN, &anVals[0]);
2520 5 : if (mappds == nullptr)
2521 : {
2522 0 : CPLError(CE_Failure, CPLE_AppDefined,
2523 : "Could not get extended template definition");
2524 0 : CSLDestroy(papszTokens);
2525 0 : return false;
2526 : }
2527 5 : if (nTokens < mappds->maplen + mappds->extlen)
2528 : {
2529 1 : CPLError(CE_Failure, CPLE_AppDefined,
2530 : "PDS_PDTN = %d (with provided elements) requires "
2531 : "%d elements in PDS_TEMPLATE_ASSEMBLED_VALUES. "
2532 : "Only %d provided",
2533 1 : nPDTN, mappds->maplen + mappds->extlen, nTokens);
2534 1 : free(mappds->ext);
2535 1 : free(mappds);
2536 1 : CSLDestroy(papszTokens);
2537 1 : return false;
2538 : }
2539 4 : else if (nTokens > mappds->maplen + mappds->extlen)
2540 : {
2541 1 : CPLError(CE_Warning, CPLE_AppDefined,
2542 : "PDS_PDTN = %d (with provided elements) requires"
2543 : "%d elements in PDS_TEMPLATE_ASSEMBLED_VALUES. "
2544 : "But %d provided. Extra elements will be ignored",
2545 1 : nPDTN, mappds->maplen + mappds->extlen, nTokens);
2546 : }
2547 :
2548 4 : WriteAssembledPDS(fp, mappds, true, papszTokens, anVals);
2549 : }
2550 :
2551 7 : free(mappds->ext);
2552 7 : free(mappds);
2553 7 : CSLDestroy(papszTokens);
2554 : }
2555 156 : PatchSectionSize(fp, nStartSection4);
2556 156 : return true;
2557 : }
2558 :
2559 : /************************************************************************/
2560 : /* CreateCopy() */
2561 : /************************************************************************/
2562 :
2563 159 : GDALDataset *GRIBDataset::CreateCopy(const char *pszFilename,
2564 : GDALDataset *poSrcDS, int /* bStrict */,
2565 : char **papszOptions,
2566 : GDALProgressFunc pfnProgress,
2567 : void *pProgressData)
2568 :
2569 : {
2570 318 : if (poSrcDS->GetRasterYSize() == 0 ||
2571 159 : poSrcDS->GetRasterXSize() > INT_MAX / poSrcDS->GetRasterXSize())
2572 : {
2573 1 : CPLError(CE_Failure, CPLE_NotSupported,
2574 : "Cannot create GRIB2 rasters with more than 2 billion pixels");
2575 1 : return nullptr;
2576 : }
2577 :
2578 158 : GDALGeoTransform gt;
2579 158 : if (poSrcDS->GetGeoTransform(gt) != CE_None)
2580 : {
2581 1 : CPLError(CE_Failure, CPLE_NotSupported,
2582 : "Source dataset must have a geotransform");
2583 1 : return nullptr;
2584 : }
2585 157 : if (gt[2] != 0.0 || gt[4] != 0.0)
2586 : {
2587 1 : CPLError(CE_Failure, CPLE_NotSupported,
2588 : "Geotransform with rotation terms not supported");
2589 1 : return nullptr;
2590 : }
2591 :
2592 312 : OGRSpatialReference oSRS;
2593 156 : oSRS.importFromWkt(poSrcDS->GetProjectionRef());
2594 156 : if (oSRS.IsProjected())
2595 : {
2596 90 : const char *pszProjection = oSRS.GetAttrValue("PROJECTION");
2597 90 : if (pszProjection == nullptr ||
2598 90 : !(EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR) ||
2599 11 : EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) ||
2600 9 : EQUAL(pszProjection, SRS_PT_MERCATOR_2SP) ||
2601 6 : EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC) ||
2602 4 : EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP) ||
2603 3 : EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) ||
2604 2 : EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA) ||
2605 1 : EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA)))
2606 : {
2607 0 : CPLError(CE_Failure, CPLE_NotSupported,
2608 : "Unsupported projection: %s",
2609 : pszProjection ? pszProjection : "");
2610 0 : return nullptr;
2611 : }
2612 : }
2613 66 : else if (!oSRS.IsGeographic())
2614 : {
2615 1 : CPLError(CE_Failure, CPLE_NotSupported,
2616 : "Unsupported or missing spatial reference system");
2617 1 : return nullptr;
2618 : }
2619 :
2620 155 : const bool bAppendSubdataset = CPLTestBool(
2621 : CSLFetchNameValueDef(papszOptions, "APPEND_SUBDATASET", "NO"));
2622 155 : VSILFILE *fp = VSIFOpenL(pszFilename, bAppendSubdataset ? "rb+" : "wb+");
2623 155 : if (fp == nullptr)
2624 : {
2625 4 : CPLError(CE_Failure, CPLE_FileIO, "Cannot create %s", pszFilename);
2626 4 : return nullptr;
2627 : }
2628 151 : VSIFSeekL(fp, 0, SEEK_END);
2629 :
2630 151 : vsi_l_offset nStartOffset = 0;
2631 151 : vsi_l_offset nTotalSizeOffset = 0;
2632 151 : int nSplitAndSwapColumn = 0;
2633 : // Note: WRITE_SUBGRIDS=YES should not be used blindly currently, as it
2634 : // does not check that the content of the DISCIPLINE and IDS are the same.
2635 : // A smarter behavior would be to break into separate messages if needed
2636 : const bool bWriteSubGrids =
2637 151 : CPLTestBool(CSLFetchNameValueDef(papszOptions, "WRITE_SUBGRIDS", "NO"));
2638 296 : for (int nBand = 1; nBand <= poSrcDS->GetRasterCount(); nBand++)
2639 : {
2640 162 : if (nBand == 1 || !bWriteSubGrids)
2641 : {
2642 : // Section 0: Indicator section
2643 162 : nStartOffset = VSIFTellL(fp);
2644 162 : VSIFWriteL("GRIB", 4, 1, fp);
2645 162 : WriteByte(fp, 0); // reserved
2646 162 : WriteByte(fp, 0); // reserved
2647 : int nDiscipline =
2648 162 : atoi(GetBandOption(papszOptions, poSrcDS, nBand, "DISCIPLINE",
2649 : "0")); // 0 = Meteorological
2650 162 : WriteByte(fp, nDiscipline); // discipline
2651 162 : WriteByte(fp, 2); // GRIB edition number
2652 162 : nTotalSizeOffset = VSIFTellL(fp);
2653 162 : WriteUInt32(fp, GRIB2MISSING_u4); // dummy file size (high 32 bits)
2654 162 : WriteUInt32(fp, GRIB2MISSING_u4); // dummy file size (low 32 bits)
2655 :
2656 : // Section 1: Identification Section
2657 162 : WriteSection1(fp, poSrcDS, nBand, papszOptions);
2658 :
2659 : // Section 2: Local use section
2660 162 : WriteUInt32(fp, 5); // section size
2661 162 : WriteByte(fp, 2); // section number
2662 :
2663 : // Section 3: Grid Definition Section
2664 162 : GRIB2Section3Writer oSection3(fp, poSrcDS);
2665 162 : if (!oSection3.Write())
2666 : {
2667 0 : VSIFCloseL(fp);
2668 0 : return nullptr;
2669 : }
2670 162 : nSplitAndSwapColumn = oSection3.SplitAndSwap();
2671 : }
2672 :
2673 : // Section 4: Product Definition Section
2674 162 : float fValOffset = 0.0f;
2675 162 : if (!WriteSection4(fp, poSrcDS, nBand, papszOptions, fValOffset))
2676 : {
2677 6 : VSIFCloseL(fp);
2678 6 : return nullptr;
2679 : }
2680 :
2681 : // Section 5, 6 and 7
2682 156 : if (!GRIB2Section567Writer(fp, poSrcDS, nBand, nSplitAndSwapColumn)
2683 156 : .Write(fValOffset, papszOptions, pfnProgress, pProgressData))
2684 : {
2685 11 : VSIFCloseL(fp);
2686 11 : return nullptr;
2687 : }
2688 :
2689 145 : if (nBand == poSrcDS->GetRasterCount() || !bWriteSubGrids)
2690 : {
2691 : // Section 8: End section
2692 145 : VSIFWriteL("7777", 4, 1, fp);
2693 :
2694 : // Patch total message size at end of section 0
2695 145 : vsi_l_offset nCurOffset = VSIFTellL(fp);
2696 145 : if (nCurOffset - nStartOffset > INT_MAX)
2697 : {
2698 0 : CPLError(CE_Failure, CPLE_NotSupported,
2699 : "GRIB message larger than 2 GB");
2700 0 : VSIFCloseL(fp);
2701 0 : return nullptr;
2702 : }
2703 145 : GUInt32 nTotalSize =
2704 145 : static_cast<GUInt32>(nCurOffset - nStartOffset);
2705 145 : VSIFSeekL(fp, nTotalSizeOffset, SEEK_SET);
2706 145 : WriteUInt32(fp, 0); // file size (high 32 bits)
2707 145 : WriteUInt32(fp, nTotalSize); // file size (low 32 bits)
2708 :
2709 145 : VSIFSeekL(fp, nCurOffset, SEEK_SET);
2710 : }
2711 :
2712 290 : if (pfnProgress &&
2713 145 : !pfnProgress(static_cast<double>(nBand) / poSrcDS->GetRasterCount(),
2714 : nullptr, pProgressData))
2715 : {
2716 0 : VSIFCloseL(fp);
2717 0 : return nullptr;
2718 : }
2719 : }
2720 :
2721 134 : VSIFCloseL(fp);
2722 :
2723 268 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
2724 134 : return Open(&oOpenInfo);
2725 : }
|
