WSL/SLF GitLab Repository

PNGIO.cc 28.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#include "PNGIO.h"
#include <meteoio/ResamplingAlgorithms2D.h>
20
#include <meteoio/Graphics.h>
21
#include <meteoio/meteolaws/Meteoconst.h>
22

23
#include <cstring>
24
#include <algorithm>
25
#include <errno.h>
26
#include <zlib.h>
27
28
29
30
31

using namespace std;

namespace mio {
/**
32
 * @page pngio PNGIO
33
 * @section template_format Format
34
35
36
37
 * This plugin write data to the Portable Network Graphics format (see https://secure.wikimedia.org/wikipedia/en/wiki/Portable_Network_Graphics).
 * No data read has been implemented, because reading an existing file would require the exact knowlege of the color gradient that has been used
 * to create it. When writing grids, various color gradients will be used depending on the parameter that the data represents. Nodata values
 * are represented by transparent pixels (transparency is acheived through a transparent color instead of a true alpha channel for size and performance).
38
 * If a grid containing no data (ie: size 0x0) is sent to the plugin, then no file will be written.
39
 * Finally, the naming scheme for meteo grids should be: YYYY-MM-DDTHH.mm_{MeteoGrids::Parameters}.png
40
41
 *
 * @section template_units Units
42
 * All units are MKSA except temperatures that are expressed in celcius.
43
44
45
 *
 * @section template_keywords Keywords
 * This plugin uses the following keywords:
46
47
 * - COORDSYS: input coordinate system (see Coords) specified in the [Output] section
 * - COORDPARAM: extra input coordinates parameters (see Coords) specified in the [Output] section
48
49
50
51
52
53
54
 * - GRID2DPATH: meteo grids directory where to read the grids; [Output] section
 * - PNG_LEGEND: plot legend on the side of the graph? (default: true)
 * - PNG_MIN_SIZE: guarantee that a 2D plot will have at least the given size
 * - PNG_MAX_SIZE: guarantee that a 2D plot will have at most the given size
 * - PNG_SCALING: scaling algorithm, either nearest or bilinear (default=bilinear)
 * - PNG_AUTOSCALE: autoscale for the color gradient? (default=true)
 * - PNG_WORLD_FILE: create world file with each file? (default=false)
55
56
 *
 * Advanced parameters (ie: don't mess up with them if you don't know what you're doing):
57
 * - PNG_INDEXED: create an indexed PNG? (default=true)
58
 * - PNG_NR_LEVELS: number of colors to use (less=smaller files, but it must be at least 5 and less than 255. default=30)
59
 * - PNG_SPEED_OPTIMIZE: optimize file creation for speed? (default=true, otherwise optimize for file size)
60
61
 *
 * The size are specified as width followed by height, with the separator being either a space, 'x' or '*'. If a minimum and a maximum size are given, the average of the smallest and largest permissible sizes will be used.
62
63
64
65
66
67
68
69
70
71
72
 * The world file is used for geolocalization and goes alongside the graphics output. By convention,
 * the file has the same name as the image file, with the third letter of the extension jammed with a w: tif->tfw, jpg->jqw.
 * The format is the following:
 * @code
 *    5.000000000000 (size of pixel in x direction)
 *    0.000000000000 (rotation term for row)
 *    0.000000000000 (rotation term for column)
 *    -5.000000000000 (size of pixel in y direction)
 *    492169.690845528910 (x coordinate of centre of upper left pixel in map units)
 *    5426523.318065105000 (y coordinate of centre of upper left pixel in map units)
 * @endcode
73
 *
74
 * @section example Example use
75
76
77
78
79
80
 * @code
 * GRID2D = PNG
 * png_legend = false
 * png_min_size = 400x400
 * png_max_size = 1366*768
 * @endcode
81
82
83
84
85
86
87
88
89
90
91
92
93
94
 *
 * @section Compilation
 * In order to compile this plugin, you need libpng and zlib. For Linux, please select both the libraries and their development files in your package manager.
 *
 * For Windows, you can find zlib at http://switch.dl.sourceforge.net/project/gnuwin32/zlib/1.2.3/zlib-1.2.3.exe
 * and libpng at http://switch.dl.sourceforge.net/project/gnuwin32/libpng/1.2.37/libpng-1.2.37-setup.exe . Once this has been installed, if you plan on using
 * Visual c++, you also need to edit the file zconf.h in the libpng installation directory and transform the line 287:
 * @code
 * #if 0           // HAVE_UNISTD_H etc etc
 * @endcode
 * should become
 * @code
 * #if 1           // HAVE_UNISTD_H etc etc
 * @endcode
95
96
97
 */

const double PNGIO::plugin_nodata = -999.; //plugin specific nodata value. It can also be read by the plugin (depending on what is appropriate)
98
99
100
const unsigned char PNGIO::channel_depth = 8;
const unsigned char PNGIO::channel_max_color = 255;
const unsigned char PNGIO::transparent_grey = channel_max_color;
101

102
PNGIO::PNGIO(const std::string& configfile)
103
       : cfg(configfile),
104
105
106
107
108
         fp(NULL), autoscale(true), has_legend(true), has_world_file(false), optimize_for_speed(true),
         indexed_png(true), nr_levels(30),
         coordout(), coordoutparam(), grid2dpath(),
         scaling("bilinear"), min_w(IOUtils::unodata), min_h(IOUtils::unodata), max_w(IOUtils::unodata), max_h(IOUtils::unodata),
         metadata_key(), metadata_text()
109
{
110
	setOptions();
111
112
}

113
PNGIO::PNGIO(const Config& cfgreader)
114
       : cfg(cfgreader),
115
116
117
118
119
         fp(NULL), autoscale(true), has_legend(true), has_world_file(false), optimize_for_speed(true),
         indexed_png(true), nr_levels(30),
         coordout(), coordoutparam(), grid2dpath(),
         scaling("bilinear"), min_w(IOUtils::unodata), min_h(IOUtils::unodata), max_w(IOUtils::unodata), max_h(IOUtils::unodata),
         metadata_key(), metadata_text()
120
{
121
122
123
	setOptions();
}

124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
PNGIO& PNGIO::operator=(const PNGIO& source) {
	if(this != &source) {
		fp = NULL;
		autoscale = source.autoscale;
		has_legend = source.has_legend;
		has_world_file = source.has_world_file;
		optimize_for_speed = source.optimize_for_speed;
		indexed_png = source.indexed_png;
		nr_levels = source.nr_levels;
		coordout = source.coordout;
		coordoutparam = source.coordoutparam;
		grid2dpath = source.grid2dpath;
		scaling = source.scaling;
		min_w = source.min_w;
		min_h = source.min_h;
		max_w = source.max_w;
		max_h = source.max_h;
		metadata_key = source.metadata_key;
		metadata_text = source.metadata_text;
	}
	return *this;
}

147
148
void PNGIO::setOptions()
{
149
	//default values have been set by the constructors
150
	cfg.getValue("COORDSYS", "Output", coordout);
151
	cfg.getValue("COORDPARAM", "Output", coordoutparam, IOUtils::nothrow);
152
	cfg.getValue("GRID2DPATH", "Output", grid2dpath);
153
	//cfg.getValue("TIME_ZONE", "Output", tz_out, IOUtils::nothrow);
154

155
156
	//get size specifications
	std::string min_size, max_size;
157
	cfg.getValue("PNG_MIN_SIZE", "Output", min_size, IOUtils::nothrow);
158
	if(!min_size.empty()) parse_size(min_size, min_w, min_h);
159
	cfg.getValue("PNG_MAX_SIZE", "Output", max_size, IOUtils::nothrow);
160
	if(!max_size.empty()) parse_size(max_size, max_w, max_h);
161

162
163
164
165
	cfg.getValue("PNG_AUTOSCALE", "Output", autoscale, IOUtils::nothrow);
	cfg.getValue("PNG_LEGEND", "Output", has_legend, IOUtils::nothrow);
	cfg.getValue("PNG_SCALING", "Output", scaling, IOUtils::nothrow);
	cfg.getValue("PNG_WORLD_FILE", "Output", has_world_file, IOUtils::nothrow);
166
167
168
169
170

	if(has_legend) { //we need to save room for the legend
		if(min_w!=IOUtils::unodata) min_w -= legend::getLegendWidth();
		if(max_w!=IOUtils::unodata) max_w -= legend::getLegendWidth();
	}
171

172
173
	cfg.getValue("PNG_INDEXED", "Output", indexed_png, IOUtils::nothrow);
	cfg.getValue("PNG_SPEED_OPTIMIZE", "Output", optimize_for_speed, IOUtils::nothrow);
174
	unsigned int tmp=IOUtils::unodata;
175
	cfg.getValue("PNG_NR_LEVELS", "Output", tmp, IOUtils::nothrow);
176
	if(tmp!=IOUtils::unodata && (tmp>255 || tmp<5)) {
177
178
		throw InvalidFormatException("PNG_NR_LEVELS must be between 5 and 255!", AT);
	}
179
	if(tmp!=IOUtils::unodata) nr_levels=static_cast<unsigned char>(tmp);
180
181
}

182
void PNGIO::parse_size(const std::string& size_spec, size_t& width, size_t& height)
183
184
{
	char rest[32] = "";
185
186
187
	if(sscanf(size_spec.c_str(), "%zu %zu%31s", &width, &height, rest) < 2)
	if(sscanf(size_spec.c_str(), "%zu*%zu%31s", &width, &height, rest) < 2)
	if(sscanf(size_spec.c_str(), "%zux%zu%31s", &width, &height, rest) < 2) {
188
		std::ostringstream ss;
189
190
191
192
193
194
		ss << "Can not parse PNGIO size specification \"" << size_spec << "\"";
		throw InvalidFormatException(ss.str(), AT);
	}
	std::string tmp(rest);
	IOUtils::trim(tmp);
	if ((tmp.length() > 0) && tmp[0] != '#' && tmp[0] != ';') {//if line holds more than one value it's invalid
195
		std::ostringstream ss;
196
197
198
199
200
		ss << "Invalid PNGIO size specification \"" << size_spec << "\"";
		throw InvalidFormatException(ss.str(), AT);
	}
}

201
double PNGIO::getScaleFactor(const size_t& grid_w, const size_t& grid_h)
202
203
204
205
206
207
208
209
210
211
212
{
	if(grid_w==0 || grid_h==0) {
		return 1.;
	}

	double min_factor = IOUtils::nodata;
	if(min_w!=IOUtils::unodata) { //min_w & min_w are read together
		const double min_w_factor = (double)min_w / (double)grid_w;
		const double min_h_factor = (double)min_h / (double)grid_h;
		min_factor = std::max(min_w_factor, min_h_factor);
	}
213

214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
	double max_factor = IOUtils::nodata;
	if(max_w!=IOUtils::unodata) { //max_w & max_h are read together
		const double max_w_factor = (double)max_w / (double)grid_w;
		const double max_h_factor = (double)max_h / (double)grid_h;
		max_factor = std::min(max_w_factor, max_h_factor);
	}

	if(min_factor==IOUtils::nodata && max_factor==IOUtils::nodata)
		return 1.; //no user given specification
	if(min_factor!=IOUtils::nodata && max_factor!=IOUtils::nodata)
		return (min_factor+max_factor)/2.; //both min & max -> average

	//only one size specification provided -> return its matching factor
	if(min_factor!=IOUtils::nodata)
		return min_factor;
	else
		return max_factor;
231
232
}

233
PNGIO::~PNGIO() throw() {
234
	if(fp!=NULL) fclose(fp); fp=NULL;
235
236
}

237
238
239
240
241
242
243
void PNGIO::read2DGrid(Grid2DObject&, const std::string&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void PNGIO::read2DGrid(Grid2DObject&, const MeteoGrids::Parameters& , const Date&)
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void PNGIO::readDEM(DEMObject& /*dem_out*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void PNGIO::readLanduse(Grid2DObject& /*landuse_out*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void PNGIO::readAssimilationData(const Date& /*date_in*/, Grid2DObject& /*da_out*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void PNGIO::readStationData(const Date&, std::vector<StationData>& /*vecStation*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void PNGIO::readMeteoData(const Date& /*dateStart*/, const Date& /*dateEnd*/,
                             std::vector< std::vector<MeteoData> >& /*vecMeteo*/,
                             const size_t&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void PNGIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& /*vecMeteo*/,
                              const std::string&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void PNGIO::readSpecialPoints(std::vector<Coords>&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

294
295
Grid2DObject PNGIO::scaleGrid(const Grid2DObject& grid_in)
{ //scale input image
296
	const double factor = getScaleFactor(grid_in.ncols, grid_in.nrows);
297
	if(scaling=="nearest")
298
		return ResamplingAlgorithms2D::NearestNeighbour(grid_in, factor);
299
	else if(scaling=="bilinear")
300
		return ResamplingAlgorithms2D::BilinearResampling(grid_in, factor);
301
	else {
302
		ostringstream ss;
303
304
305
		ss << "Grid scaling algorithm \"" << scaling << "\" unknown";
		throw UnknownValueException(ss.str(), AT);
	}
306
}
307

308
void PNGIO::setFile(const std::string& filename, png_structp& png_ptr, png_infop& info_ptr, const size_t &width, const size_t &height)
309
{
310
311
312
313
	// Open file for writing (binary mode)
	if (!IOUtils::validFileName(filename)) {
		throw InvalidFileNameException(filename, AT);
	}
314
	errno=0;
315
316
	fp = fopen(filename.c_str(), "wb");
	if (fp == NULL) {
317
		ostringstream ss;
318
319
		ss << "Error openning file \"" << filename << "\", possible reason: " << strerror(errno);
		throw FileAccessException(ss.str(), AT);
320
321
322
323
324
	}

	// Initialize write structure
	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
	if (png_ptr == NULL) {
325
		fclose(fp); fp=NULL;
326
327
328
329
330
331
		throw IOException("Could not allocate write structure", AT);
	}

	// Initialize info structure
	info_ptr = png_create_info_struct(png_ptr);
	if (info_ptr == NULL) {
332
		fclose(fp); fp=NULL;
333
		png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
334
		free(png_ptr);
335
336
337
338
		throw IOException("Could not allocate info structure", AT);
	}

	// Setup Exception handling
339
#ifdef WIN32
340
341
	#pragma warning(disable:4611) //the setjmp of libpng has been set up so that it can safely be called from c++
#endif
342
	if (setjmp(png_jmpbuf(png_ptr))) {
343
		closePNG(png_ptr, info_ptr, NULL);
344
		throw IOException("Error during png creation. Can not set jump pointer (I have no clue what it means too!)", AT);
345
346
347
348
	}

	png_init_io(png_ptr, fp);

349
	if(optimize_for_speed) png_set_compression_level(png_ptr, Z_BEST_SPEED);
350
351
	else png_set_compression_level(png_ptr, Z_BEST_COMPRESSION);

352
353
	png_set_filter(png_ptr, PNG_FILTER_TYPE_BASE, PNG_FILTER_SUB|PNG_FILTER_UP); //any other filter is costly and brings close to nothing...
	if(indexed_png) png_set_compression_strategy(png_ptr, Z_RLE); //Z_DEFAULT_STRATEGY, Z_FILTERED, Z_HUFFMAN_ONLY, Z_RLE
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371

	// Write header (8 bit colour depth). Full alpha channel with PNG_COLOR_TYPE_RGB_ALPHA
	if(indexed_png) {
		png_set_IHDR(png_ptr, info_ptr, width, height,
			channel_depth, PNG_COLOR_TYPE_PALETTE, PNG_INTERLACE_NONE,
			PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
		//set transparent color (ie: cheap transparency: leads to smaller files and shorter run times)
		png_byte trans = 0; //by convention, the gradient define it as color 0
		png_set_tRNS(png_ptr, info_ptr, &trans, 1, 0);
	} else {
		png_set_IHDR(png_ptr, info_ptr, width, height,
			channel_depth, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
			PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
		//set transparent color (ie: cheap transparency: leads to smaller files and shorter run times)
		png_color_16 trans_rgb_value = {transparent_grey, transparent_grey, transparent_grey, transparent_grey, transparent_grey};
		png_set_tRNS(png_ptr, info_ptr, 0, 0, &trans_rgb_value);
	}

372
373
374
	//set background color to help applications show the picture when no background is present
	png_color_16 background = {channel_max_color, channel_max_color, channel_max_color, channel_max_color, channel_max_color};
	png_set_background(png_ptr, &background, PNG_BACKGROUND_GAMMA_SCREEN, true, 1.0);
375
376
}

377
size_t PNGIO::setLegend(const size_t &ncols, const size_t &nrows, const double &min, const double &max, Array2D<double> &legend_array)
378
{
379
	if(has_legend) {
380
		const legend leg(static_cast<unsigned int>(nrows), min, max);
381
		legend_array = leg.getLegend();
382
		size_t nx, ny;
383
		legend_array.size(nx,ny);
384
385
386
		return (ncols+nx);
	} else {
		return ncols;
387
	}
388
389
}

390
void PNGIO::writeDataSection(const Grid2DObject& grid, const Array2D<double>& legend_array, const Gradient& gradient, const size_t& full_width, const png_structp& png_ptr, png_infop& info_ptr)
391
{
392
393
	const size_t ncols = grid.ncols;
	const size_t nrows = grid.nrows;
394

395
	// Allocate memory for one row (3 bytes per pixel - RGB)
396
397
398
399
	unsigned char channels;
	if(indexed_png)
		channels = 1;
	else
400
		channels = 3; //4 for rgba
401

402
403
	png_bytep row = (png_bytep)calloc(channels*sizeof(png_byte), full_width);
	if(row==NULL) {
404
405
406
		fclose(fp); fp=NULL;
		png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
		free(png_ptr);
407
408
		throw IOException("Can not allocate row memory in PNGIO!", AT);
	}
409
410

	// Write image data
411
	if(indexed_png) {
412
		for(size_t y=(nrows-1) ; y-- > 0; ) {
413
			size_t x=0;
414
			for(; x<ncols ; x++) {
415
				const size_t i=x*channels;
416
				unsigned char index;
417
				gradient.getColor(grid(x,y), index);
418
				row[i]=static_cast<png_byte>(index);
419
420
			}
			for(; x<full_width; x++) {
421
				const size_t i=x*channels;
422
				unsigned char index;
423
				gradient.getColor(legend_array(x-ncols,y), index);
424
				row[i]=static_cast<png_byte>(index);
425
			}
426
			png_write_row(png_ptr, row);
427
		}
428
	} else {
429
		for(size_t y=(nrows-1) ; y -- > 0; ) {
430
			size_t x=0;
431
			for(; x<ncols ; x++) {
432
				const size_t i=x*channels;
433
434
435
436
				unsigned char r,g,b;
				bool a;
				gradient.getColor(grid(x,y), r,g,b,a);
				if(a==true) {
437
					row[i]=static_cast<png_byte>(transparent_grey); row[i+1]=static_cast<png_byte>(transparent_grey); row[i+2]=static_cast<png_byte>(transparent_grey);
438
				} else {
439
					row[i]=static_cast<png_byte>(r); row[i+1]=static_cast<png_byte>(g); row[i+2]=static_cast<png_byte>(b);
440
441
442
				}
			}
			for(; x<full_width; x++) {
443
				const size_t i=x*channels;
444
445
446
447
				unsigned char r,g,b;
				bool a;
				gradient.getColor(legend_array(x-ncols,y), r,g,b,a);
				if(a==true) {
448
					row[i]=static_cast<png_byte>(transparent_grey); row[i+1]=static_cast<png_byte>(transparent_grey); row[i+2]=static_cast<png_byte>(transparent_grey);
449
				} else {
450
					row[i]=static_cast<png_byte>(r); row[i+1]=static_cast<png_byte>(g); row[i+2]=static_cast<png_byte>(b);
451
				}
452
			}
453
			png_write_row(png_ptr, row);
454
455
		}
	}
456

457
458
	png_write_flush(png_ptr);
	png_free(png_ptr, row);
459
460
}

461
void PNGIO::setPalette(const Gradient &gradient, png_structp& png_ptr, png_infop& info_ptr, png_color *palette)
462
{
463
464
465
	std::vector<unsigned char> pal;
	size_t nr_colors;
	gradient.getPalette(pal, nr_colors);
466
	palette = (png_color*)calloc(sizeof (png_color), nr_colors); //ie: three png_bytes, each being an unsigned char
467
	for(size_t ii=0; ii<nr_colors; ii++) {
468
469
470
471
		const size_t interlace = ii*3; //colors from Gradient interlaced
		palette[ii].red = static_cast<png_byte>(pal[interlace]);
		palette[ii].green = static_cast<png_byte>(pal[interlace+1]);
		palette[ii].blue = static_cast<png_byte>(pal[interlace+2]);
472
	}
473
	png_set_PLTE(png_ptr, info_ptr, palette, static_cast<int>(nr_colors));
474
475
}

476
void PNGIO::closePNG(png_structp& png_ptr, png_infop& info_ptr, png_color *palette)
477
{
478
	png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
479
	if(indexed_png && palette!=NULL) free(palette);
480
	png_destroy_write_struct(&png_ptr, &info_ptr);
481
	fclose(fp); fp=NULL;
482
483
	free(info_ptr);
	free(png_ptr);
484
485
}

486
void PNGIO::write2DGrid(const Grid2DObject& grid_in, const std::string& filename)
487
{
488
	const string full_name = grid2dpath+"/"+filename;
489
	fp=NULL;
490
	png_color *palette=NULL;
491
492
493
494
495
	png_structp png_ptr=NULL;
	png_infop info_ptr=NULL;

	//scale input image
	const Grid2DObject grid = scaleGrid(grid_in);
496
	const size_t ncols = grid_in.ncols, nrows = grid_in.nrows;
497
498
	if(ncols==0 || nrows==0) return;

499
500
501
502
	const double min = grid.grid2D.getMin();
	const double max = grid.grid2D.getMax();

	Gradient gradient(Gradient::heat, min, max, autoscale);
503
	if(indexed_png) gradient.setNrOfLevels(nr_levels);
504
505

	Array2D<double> legend_array; //it will remain empty if there is no legend
506
	const size_t full_width = setLegend(ncols, nrows, min, max, legend_array);
507

508
	setFile(full_name, png_ptr, info_ptr, full_width, nrows);
509
	if(indexed_png) setPalette(gradient, png_ptr, info_ptr, palette);
510
	if(has_world_file) writeWorldFile(grid, full_name);
511
512
513
514
515
516

	createMetadata(grid);
	metadata_key.push_back("Title"); //adding generic title
	metadata_text.push_back("Unknown Gridded data");
	writeMetadata(png_ptr, info_ptr);

517
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
518
	png_write_end(png_ptr, NULL);
519

520
	closePNG(png_ptr, info_ptr, palette);
521
522
}

523
void PNGIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameters& parameter, const Date& date)
524
{
525
526
	std::string filename;
	if(parameter==MeteoGrids::DEM || parameter==MeteoGrids::SLOPE || parameter==MeteoGrids::AZI)
527
		filename = grid2dpath + "/" + MeteoGrids::getParameterName(parameter) + ".png";
528
529
530
531
532
	else {
		std::string date_str = date.toString(Date::ISO);
		std::replace( date_str.begin(), date_str.end(), ':', '.');
		filename = grid2dpath + "/" + date_str + "_" + MeteoGrids::getParameterName(parameter) + ".png";
	}
533

534
	fp=NULL;
535
	png_color *palette=NULL;
536
537
538
539
540
	png_structp png_ptr=NULL;
	png_infop info_ptr=NULL;

	//scale input image
	Grid2DObject grid = scaleGrid(grid_in);
541
	const size_t ncols = grid_in.ncols, nrows = grid_in.nrows;
542
543
	if(ncols==0 || nrows==0) return;

544
545
546
547
548
	double min = grid.grid2D.getMin();
	double max = grid.grid2D.getMax();

	Gradient gradient;
	if(parameter==MeteoGrids::DEM) {
549
550
		if(!autoscale) {
			min = 0.; //we want a 3000 snow line with a full scale legend
551
552
			max = 3500.;
			gradient.set(Gradient::terrain, min, max, autoscale); //max used as snow line reference
553
554
		} else
			gradient.set(Gradient::terrain, min, max, autoscale);
555
556
557
	} else if(parameter==MeteoGrids::SLOPE) {
		gradient.set(Gradient::slope, min, max, autoscale);
	} else if(parameter==MeteoGrids::AZI) {
558
559
560
561
		if(!autoscale) {
			min = 0.;
			max = 360.;
		}
562
		gradient.set(Gradient::azi, min, max, autoscale);
563
564
565
566
567
568
	} else if(parameter==MeteoGrids::DW) {
		if(!autoscale) {
			min = 0.;
			max = 360.;
		}
		gradient.set(Gradient::azi, min, max, autoscale);
569
	} else if(parameter==MeteoGrids::HS) {
570
		if(!autoscale) {
571
			min = 0.; max = 2.5;
572
573
		}
		gradient.set(Gradient::blue, min, max, autoscale);
574
575
	} else if(parameter==MeteoGrids::TA) {
		grid.grid2D -= Cst::t_water_freezing_pt; //convert to celsius
576
		if(!autoscale) {
577
			min = -15.; max = 15.;
578
579
580
581
		} else {
			min -= Cst::t_water_freezing_pt;
			max -= Cst::t_water_freezing_pt;
		}
582
		gradient.set(Gradient::heat, min, max, autoscale);
583
584
585
	} else if(parameter==MeteoGrids::TSS) {
		grid.grid2D -= Cst::t_water_freezing_pt; //convert to celsius
		if(!autoscale) {
586
			min = -20.; max = 5.;
587
588
589
590
591
		} else {
			min -= Cst::t_water_freezing_pt;
			max -= Cst::t_water_freezing_pt;
		}
		gradient.set(Gradient::freeze, min, max, autoscale);
592
	} else if(parameter==MeteoGrids::RH) {
593
594
595
		if(!autoscale) {
			min = 0.; max = 1.;
		}
596
		gradient.set(Gradient::bg_isomorphic, min, max, autoscale);
597
598
599
600
601
602
603
604
605
606
607
	} else if(parameter==MeteoGrids::P) {
		if(!autoscale) {
			//lowest and highest pressures ever recorded on Earth: 87000 and 108570
			min = 87000.; max = 115650.; //centered around 1 atm
			gradient.set(Gradient::bluewhitered, min, max, autoscale);
		} else {
			const double delta1 = fabs(Cst::std_press-min);
			const double delta2 = fabs(max - Cst::std_press);
			const double delta = (delta1>delta2)?delta1:delta2;
			gradient.set(Gradient::bluewhitered, Cst::std_press-delta, Cst::std_press+delta, autoscale);
		}
608
609
610
611
	} else if(parameter==MeteoGrids::ALB) {
		if(!autoscale) {
			min = 0.; max = 1.;
		}
612
		gradient.set(Gradient::blktowhite, min, max, autoscale);
613
614
	} else if(parameter==MeteoGrids::ISWR) {
		if(!autoscale) {
615
			min = 0.; max = 800.;
616
617
618
619
620
621
622
		}
		gradient.set(Gradient::heat, min, max, autoscale);
	} else if(parameter==MeteoGrids::ILWR) {
		if(!autoscale) {
			min = 200.; max = 500.;
		}
		gradient.set(Gradient::heat, min, max, autoscale);
623
	} else if(parameter==MeteoGrids::SWE) {
624
		if(!autoscale) {
625
			min = 0.; max = 250.;
626
		}
627
		gradient.set(Gradient::blue_pink, min, max, autoscale);
628
629
630
	} else {
		gradient.set(Gradient::heat, min, max, autoscale);
	}
631
	gradient.setNrOfLevels(nr_levels);
632
633

	Array2D<double> legend_array; //it will remain empty if there is no legend
634
	const size_t full_width = setLegend(ncols, nrows, min, max, legend_array);
635

636
	setFile(filename, png_ptr, info_ptr, full_width, nrows);
637
	if(indexed_png) setPalette(gradient, png_ptr, info_ptr, palette);
638
	if(has_world_file) writeWorldFile(grid, filename);
639
640
641
642

	createMetadata(grid);
	metadata_key.push_back("Title"); //adding title
	metadata_text.push_back( MeteoGrids::getParameterName(parameter)+" on "+date.toString(Date::ISO) );
643
644
	metadata_key.push_back("Simulation Date");
	metadata_text.push_back( date.toString(Date::ISO) );
645
646
	metadata_key.push_back("Simulation Parameter");
	metadata_text.push_back( MeteoGrids::getParameterName(parameter) );
647
648
	writeMetadata(png_ptr, info_ptr);

649
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
650
	png_write_end(png_ptr, NULL);
651

652
	closePNG(png_ptr, info_ptr, palette);
653
654
655
656
657
658
}

void PNGIO::writeWorldFile(const Grid2DObject& grid_in, const std::string& filename)
{
	const string world_file = IOUtils::removeExtension(filename)+".pnw";
	const double cellsize = grid_in.cellsize;
659
	Coords world_ref = grid_in.llcorner;
660
	world_ref.setProj(coordout, coordoutparam);
661
	world_ref.moveByXY(.5*cellsize, (double(grid_in.nrows)+.5)*cellsize); //moving to center of upper left cell
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681

	std::ofstream fout;
	fout.open(world_file.c_str());
	if (fout.fail()) {
		throw FileAccessException(world_file, AT);
	}

	try {
		fout << std::setprecision(12) << cellsize << "\n";
		fout << "0.000000000000\n";
		fout << "0.000000000000\n";
		fout << std::setprecision(12) << -cellsize << "\n";
		fout << std::setprecision(12) << world_ref.getEasting() << "\n";
		fout << std::setprecision(12) << world_ref.getNorthing() << "\n";
	} catch(...) {
		fout.close();
		throw FileAccessException("Failed when writing to PNG world file \""+world_file+"\"", AT);
	}

	fout.close();
682
683
}

684
void PNGIO::createMetadata(const Grid2DObject& grid)
685
{
686
687
	const double lat = grid.llcorner.getLat();
	const double lon = grid.llcorner.getLon();
688
	ostringstream ss;
689

690
691
692
	metadata_key.clear();
	metadata_text.clear();

693
694
695
696
	metadata_key.push_back("Creation Time");
	Date cr_date;
	cr_date.setFromSys();
	metadata_text.push_back( cr_date.toString(Date::ISO) );
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
	metadata_key.push_back("Author");
	metadata_text.push_back(IOUtils::getLogName());
	metadata_key.push_back("Software");
	metadata_text.push_back("MeteoIO "+getLibVersion());
	metadata_key.push_back("Position");
	metadata_text.push_back("llcorner");
	metadata_key.push_back("Cellsize");
	ss.str(""); ss << fixed << setprecision(2) << grid.cellsize;
	metadata_text.push_back(ss.str());
	metadata_key.push_back("Latitude");
	ss.str(""); ss << fixed << setprecision(6) << lat;
	metadata_text.push_back(ss.str());
	metadata_key.push_back("Longitude");
	ss.str(""); ss << fixed << setprecision(6) << lon;
	metadata_text.push_back(ss.str());
712

713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
	if(lat<0.) {
		metadata_key.push_back("LatitudeRef");
		metadata_text.push_back("S");
		metadata_key.push_back("GPSLatitude");
		metadata_text.push_back(decimal_to_dms(-lat));
	} else {
		metadata_key.push_back("LatitudeRef");
		metadata_text.push_back("N");
		metadata_key.push_back("GPSLatitude");
		metadata_text.push_back(decimal_to_dms(lat));
	}
	if(lon<0.) {
		metadata_key.push_back("LongitudeRef");
		metadata_text.push_back("W");
		metadata_key.push_back("GPSLongitude");
		metadata_text.push_back(decimal_to_dms(-lon));
	} else {
		metadata_key.push_back("LongitudeRef");
		metadata_text.push_back("E");
		metadata_key.push_back("GPSLongitude");
		metadata_text.push_back(decimal_to_dms(lon));
	}
}

void PNGIO::writeMetadata(png_structp &png_ptr, png_infop &info_ptr)
{
739
	const size_t max_len = 79; //according to the official specs' recommendation
740
	const size_t nr = metadata_key.size();
741
742
743
	png_text *info_text = (png_text *)calloc(sizeof(png_text), nr);
	char **key = (char**)calloc(sizeof(char)*max_len, nr);
	char **text = (char**)calloc(sizeof(char)*max_len, nr);
744
745

	for(size_t ii=0; ii<nr; ii++) {
746
747
748
749
		key[ii] = (char *)calloc(sizeof(char), max_len);
		text[ii] = (char *)calloc(sizeof(char), max_len);
		strncpy(key[ii], metadata_key[ii].c_str(), max_len);
		strncpy(text[ii], metadata_text[ii].c_str(), max_len);
750
751
752
753
		info_text[ii].key = key[ii];
		info_text[ii].text = text[ii];
		info_text[ii].compression = PNG_TEXT_COMPRESSION_NONE;
	}
754

755
	png_set_text(png_ptr, info_ptr, info_text, static_cast<int>(nr));
756
757
	png_write_info(png_ptr, info_ptr);

758
759
760
761
762
763
764
765
766
767
768
769
770
771
	free(info_text);
	for(size_t ii=0; ii<nr; ii++) {
		free(key[ii]);
		free(text[ii]);
	}
	free(key);
	free(text);
}

std::string PNGIO::decimal_to_dms(const double& decimal) {
	const int d = static_cast<int>( floor(decimal) );
	const double m = floor( ((decimal - (double)d)*60.)*100. ) / 100.;
	const double s = 3600.*(decimal - (double)d) - 60.*m;

772
	std::ostringstream dms;
773
774
	dms << d << "/1 " << static_cast<int>(m*100) << "/100 " << fixed << setprecision(6) << s << "/1";
	return dms.str();
775
776
777
}

} //namespace