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