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