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