WSL/SLF GitLab Repository

PNGIO.cc 29 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 pngio_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 pngio_units Units
42
 * All units are MKSA except temperatures that are expressed in celcius.
43
 *
44
 * @section pngio_keywords Keywords
45
 * 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 pngio_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
 * @section pngio_compilation
83
84
 * 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.
 *
85
86
87
88
89
 * For Mac, you can either install using <A href="http://www.finkproject.org/">Fink</A> or directly from <a href="http://www.libpng.org">source</a>. 
 * In this case, please install <a href="http://zlib.net/">zlib</a> at the prefix of your choice (for example with the <i>"--prefix=/usr/local"</i> options 
 * to its configure script, build with <i>"make"</i> and install with <i>"sudo make install"</i>). Then provide the libpng configure script with the 
 * <i>"--enable-shared --prefix=/usr/local"</i> options (and build with <i>"make"</i> followed by <i>"sudo make install"</i>).
 * 
90
91
92
93
94
95
96
97
98
99
 * 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
100
101
102
 */

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

107
PNGIO::PNGIO(const std::string& configfile)
108
       : cfg(configfile),
109
110
111
112
113
         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()
114
{
115
	setOptions();
116
117
}

118
PNGIO::PNGIO(const Config& cfgreader)
119
       : cfg(cfgreader),
120
121
122
123
124
         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()
125
{
126
127
128
	setOptions();
}

129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
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;
}

152
153
void PNGIO::setOptions()
{
154
	//default values have been set by the constructors
155
	cfg.getValue("COORDSYS", "Output", coordout);
156
	cfg.getValue("COORDPARAM", "Output", coordoutparam, IOUtils::nothrow);
157
	cfg.getValue("GRID2DPATH", "Output", grid2dpath);
158
	//cfg.getValue("TIME_ZONE", "Output", tz_out, IOUtils::nothrow);
159

160
161
	//get size specifications
	std::string min_size, max_size;
162
	cfg.getValue("PNG_MIN_SIZE", "Output", min_size, IOUtils::nothrow);
163
	if(!min_size.empty()) parse_size(min_size, min_w, min_h);
164
	cfg.getValue("PNG_MAX_SIZE", "Output", max_size, IOUtils::nothrow);
165
	if(!max_size.empty()) parse_size(max_size, max_w, max_h);
166

167
168
169
170
	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);
171
172
173
174
175

	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();
	}
176

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

187
void PNGIO::parse_size(const std::string& size_spec, size_t& width, size_t& height)
188
189
{
	char rest[32] = "";
190
191
192
193
	unsigned int w,h;
	if(sscanf(size_spec.c_str(), "%u %u%31s", &w, &h, rest) < 2)
	if(sscanf(size_spec.c_str(), "%u*%u%31s", &w, &h, rest) < 2)
	if(sscanf(size_spec.c_str(), "%ux%u%31s", &w, &h, rest) < 2) {
194
		std::ostringstream ss;
195
196
197
		ss << "Can not parse PNGIO size specification \"" << size_spec << "\"";
		throw InvalidFormatException(ss.str(), AT);
	}
198
199
200
201

	width = static_cast<size_t>(w);
	height = static_cast<size_t>(h);

202
203
204
	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
205
		std::ostringstream ss;
206
207
208
209
210
		ss << "Invalid PNGIO size specification \"" << size_spec << "\"";
		throw InvalidFormatException(ss.str(), AT);
	}
}

211
double PNGIO::getScaleFactor(const size_t& grid_w, const size_t& grid_h)
212
213
214
215
216
217
218
219
220
221
222
{
	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);
	}
223

224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
	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;
241
242
}

243
PNGIO::~PNGIO() throw() {
244
	if(fp!=NULL) fclose(fp); fp=NULL;
245
246
}

247
248
249
250
251
252
253
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&)
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
294
295
296
297
{
	//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);
}

298
void PNGIO::readPOI(std::vector<Coords>&)
299
300
301
302
303
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

304
305
Grid2DObject PNGIO::scaleGrid(const Grid2DObject& grid_in)
{ //scale input image
306
	const double factor = getScaleFactor(grid_in.getNx(), grid_in.getNy());
307
	if(scaling=="nearest")
308
		return ResamplingAlgorithms2D::NearestNeighbour(grid_in, factor);
309
	else if(scaling=="bilinear")
310
		return ResamplingAlgorithms2D::BilinearResampling(grid_in, factor);
311
	else {
312
		ostringstream ss;
313
314
315
		ss << "Grid scaling algorithm \"" << scaling << "\" unknown";
		throw UnknownValueException(ss.str(), AT);
	}
316
}
317

318
void PNGIO::setFile(const std::string& filename, png_structp& png_ptr, png_infop& info_ptr, const size_t &width, const size_t &height)
319
{
320
321
322
323
	// Open file for writing (binary mode)
	if (!IOUtils::validFileName(filename)) {
		throw InvalidFileNameException(filename, AT);
	}
324
	errno=0;
325
326
	fp = fopen(filename.c_str(), "wb");
	if (fp == NULL) {
327
		ostringstream ss;
328
		ss << "Error opening file \"" << filename << "\", possible reason: " << strerror(errno);
329
		throw FileAccessException(ss.str(), AT);
330
331
332
333
334
	}

	// Initialize write structure
	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
	if (png_ptr == NULL) {
335
		fclose(fp); fp=NULL;
336
		throw IOException("Could not allocate write structure. Are you running with the same version of libpng that you linked to?", AT);
337
338
339
340
341
	}

	// Initialize info structure
	info_ptr = png_create_info_struct(png_ptr);
	if (info_ptr == NULL) {
342
		fclose(fp); fp=NULL;
343
		png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
344
		free(png_ptr);
345
346
347
348
		throw IOException("Could not allocate info structure", AT);
	}

	// Setup Exception handling
349
#ifdef _MSC_VER
350
351
	#pragma warning(disable:4611) //the setjmp of libpng has been set up so that it can safely be called from c++
#endif
352
	if (setjmp(png_jmpbuf(png_ptr))) {
353
		closePNG(png_ptr, info_ptr, NULL);
354
		throw IOException("Error during png creation. Can not set jump pointer (I have no clue what it means too!)", AT);
355
356
357
358
	}

	png_init_io(png_ptr, fp);

359
	if(optimize_for_speed) png_set_compression_level(png_ptr, Z_BEST_SPEED);
360
361
	else png_set_compression_level(png_ptr, Z_BEST_COMPRESSION);

362
363
	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
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381

	// 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);
	}

382
383
384
	//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);
385
386
}

387
size_t PNGIO::setLegend(const size_t &ncols, const size_t &nrows, const double &min, const double &max, Array2D<double> &legend_array)
388
{
389
	if(has_legend) {
390
		const legend leg(static_cast<unsigned int>(nrows), min, max);
391
		legend_array = leg.getLegend();
392
		const size_t nx = legend_array.getNx();
393
394
395
		return (ncols+nx);
	} else {
		return ncols;
396
	}
397
398
}

399
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)
400
{
401
402
	const size_t ncols = grid.getNx();
	const size_t nrows = grid.getNy();
403

404
	// Allocate memory for one row (3 bytes per pixel - RGB)
405
	const unsigned char channels = (indexed_png)? 1 : 3; //4 for rgba
406
407
	png_bytep row = (png_bytep)calloc(channels*sizeof(png_byte), full_width);
	if(row==NULL) {
408
409
410
		fclose(fp); fp=NULL;
		png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
		free(png_ptr);
411
		throw IOException("Can not allocate row memory for PNG!", AT);
412
	}
413
414

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

461
462
	png_write_flush(png_ptr);
	png_free(png_ptr, row);
463
464
}

465
void PNGIO::setPalette(const Gradient &gradient, png_structp& png_ptr, png_infop& info_ptr, png_color *palette)
466
{
467
468
469
	std::vector<unsigned char> pal;
	size_t nr_colors;
	gradient.getPalette(pal, nr_colors);
470
	palette = (png_color*)calloc(sizeof (png_color), nr_colors); //ie: three png_bytes, each being an unsigned char
471
	for(size_t ii=0; ii<nr_colors; ii++) {
472
473
474
475
		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]);
476
	}
477
	png_set_PLTE(png_ptr, info_ptr, palette, static_cast<int>(nr_colors));
478
479
}

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

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

	//scale input image
	const Grid2DObject grid = scaleGrid(grid_in);
500
	const size_t ncols = grid.getNx(), nrows = grid.getNy();
501
502
	if(ncols==0 || nrows==0) return;

503
504
505
506
	const double min = grid.grid2D.getMin();
	const double max = grid.grid2D.getMax();

	Gradient gradient(Gradient::heat, min, max, autoscale);
507
	if(indexed_png) gradient.setNrOfLevels(nr_levels);
508
509

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

512
	setFile(full_name, png_ptr, info_ptr, full_width, nrows);
513
	if(indexed_png) setPalette(gradient, png_ptr, info_ptr, palette);
514
	if(has_world_file) writeWorldFile(grid, full_name);
515
516
517
518
519
520

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

521
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
522
	png_write_end(png_ptr, NULL);
523

524
	closePNG(png_ptr, info_ptr, palette);
525
526
}

527
void PNGIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameters& parameter, const Date& date)
528
{
529
530
	std::string filename;
	if(parameter==MeteoGrids::DEM || parameter==MeteoGrids::SLOPE || parameter==MeteoGrids::AZI)
531
		filename = grid2dpath + "/" + MeteoGrids::getParameterName(parameter) + ".png";
532
533
534
535
536
	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";
	}
537

538
	fp=NULL;
539
	png_color *palette=NULL;
540
541
542
543
544
	png_structp png_ptr=NULL;
	png_infop info_ptr=NULL;

	//scale input image
	Grid2DObject grid = scaleGrid(grid_in);
545
	const size_t ncols = grid.getNx(), nrows = grid.getNy();
546
547
	if(ncols==0 || nrows==0) return;

548
549
550
551
552
	double min = grid.grid2D.getMin();
	double max = grid.grid2D.getMax();

	Gradient gradient;
	if(parameter==MeteoGrids::DEM) {
553
554
		if(!autoscale) {
			min = 0.; //we want a 3000 snow line with a full scale legend
555
556
			max = 3500.;
			gradient.set(Gradient::terrain, min, max, autoscale); //max used as snow line reference
557
558
		} else
			gradient.set(Gradient::terrain, min, max, autoscale);
559
560
	} else if(parameter==MeteoGrids::SLOPE) {
		gradient.set(Gradient::slope, min, max, autoscale);
561
562
563
564
565
	} else if(parameter==MeteoGrids::SHADE) {
		if(!autoscale) {
			min = 0.; max = 1.;
		}
		gradient.set(Gradient::blktowhite, min, max, autoscale);
566
	} else if(parameter==MeteoGrids::AZI) {
567
568
569
570
		if(!autoscale) {
			min = 0.;
			max = 360.;
		}
571
		gradient.set(Gradient::azi, min, max, autoscale);
572
573
574
575
576
577
	} else if(parameter==MeteoGrids::DW) {
		if(!autoscale) {
			min = 0.;
			max = 360.;
		}
		gradient.set(Gradient::azi, min, max, autoscale);
578
	} else if(parameter==MeteoGrids::HS) {
579
		if(!autoscale) {
580
			min = 0.; max = 2.5;
581
582
		}
		gradient.set(Gradient::blue, min, max, autoscale);
583
584
	} else if(parameter==MeteoGrids::TA) {
		grid.grid2D -= Cst::t_water_freezing_pt; //convert to celsius
585
		if(!autoscale) {
586
			min = -15.; max = 15.;
587
588
589
590
		} else {
			min -= Cst::t_water_freezing_pt;
			max -= Cst::t_water_freezing_pt;
		}
591
		gradient.set(Gradient::heat, min, max, autoscale);
592
593
594
	} else if(parameter==MeteoGrids::TSS) {
		grid.grid2D -= Cst::t_water_freezing_pt; //convert to celsius
		if(!autoscale) {
595
			min = -20.; max = 5.;
596
597
598
599
600
		} else {
			min -= Cst::t_water_freezing_pt;
			max -= Cst::t_water_freezing_pt;
		}
		gradient.set(Gradient::freeze, min, max, autoscale);
601
	} else if(parameter==MeteoGrids::RH) {
602
603
604
		if(!autoscale) {
			min = 0.; max = 1.;
		}
605
		gradient.set(Gradient::bg_isomorphic, min, max, autoscale);
606
607
608
609
610
611
612
613
614
615
616
	} 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);
		}
617
618
619
620
	} else if(parameter==MeteoGrids::ALB) {
		if(!autoscale) {
			min = 0.; max = 1.;
		}
621
		gradient.set(Gradient::blktowhite, min, max, autoscale);
622
	} else if(parameter==MeteoGrids::TAU_CLD) {
623
624
625
626
		if(!autoscale) {
			min = 0.; max = 1.;
		}
		gradient.set(Gradient::blktowhite, min, max, autoscale);
627
	} else if(parameter==MeteoGrids::ISWR) {
628
		if(!autoscale) {
629
			min = 0.; max = 800.;
630
631
632
633
634
635
636
		}
		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);
637
	} else if(parameter==MeteoGrids::SWE) {
638
		if(!autoscale) {
639
			min = 0.; max = 250.;
640
		}
641
		gradient.set(Gradient::blue_pink, min, max, autoscale);
642
643
644
	} else {
		gradient.set(Gradient::heat, min, max, autoscale);
	}
645
	gradient.setNrOfLevels(nr_levels);
646
647

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

650
	setFile(filename, png_ptr, info_ptr, full_width, nrows);
651
	if(indexed_png) setPalette(gradient, png_ptr, info_ptr, palette);
652
	if(has_world_file) writeWorldFile(grid, filename);
653
654
655
656

	createMetadata(grid);
	metadata_key.push_back("Title"); //adding title
	metadata_text.push_back( MeteoGrids::getParameterName(parameter)+" on "+date.toString(Date::ISO) );
657
658
	metadata_key.push_back("Simulation Date");
	metadata_text.push_back( date.toString(Date::ISO) );
659
660
	metadata_key.push_back("Simulation Parameter");
	metadata_text.push_back( MeteoGrids::getParameterName(parameter) );
661
662
	writeMetadata(png_ptr, info_ptr);

663
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
664
	png_write_end(png_ptr, NULL);
665

666
	closePNG(png_ptr, info_ptr, palette);
667
668
669
670
671
672
}

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;
673
	Coords world_ref = grid_in.llcorner;
674
	world_ref.setProj(coordout, coordoutparam);
675
	world_ref.moveByXY(.5*cellsize, (double(grid_in.getNy())+.5)*cellsize); //moving to center of upper left cell
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695

	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();
696
697
}

698
void PNGIO::createMetadata(const Grid2DObject& grid)
699
{
700
701
	const double lat = grid.llcorner.getLat();
	const double lon = grid.llcorner.getLon();
702
	ostringstream ss;
703

704
705
706
	metadata_key.clear();
	metadata_text.clear();

707
708
709
710
	metadata_key.push_back("Creation Time");
	Date cr_date;
	cr_date.setFromSys();
	metadata_text.push_back( cr_date.toString(Date::ISO) );
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
	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());
726

727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
	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)
{
753
	const size_t max_len = 79; //according to the official specs' recommendation
754
	const size_t nr = metadata_key.size();
755
756
757
	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);
758
759

	for(size_t ii=0; ii<nr; ii++) {
760
761
762
763
		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);
764
765
766
767
		info_text[ii].key = key[ii];
		info_text[ii].text = text[ii];
		info_text[ii].compression = PNG_TEXT_COMPRESSION_NONE;
	}
768

769
	png_set_text(png_ptr, info_ptr, info_text, static_cast<int>(nr));
770
771
	png_write_info(png_ptr, info_ptr);

772
773
774
775
776
777
778
779
780
781
782
783
784
785
	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;

786
	std::ostringstream dms;
787
788
	dms << d << "/1 " << static_cast<int>(m*100) << "/100 " << fixed << setprecision(6) << s << "/1";
	return dms.str();
789
790
791
}

} //namespace