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