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