WSL/SLF GitLab Repository

PNGIO.cc 28.7 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
PNGIO& PNGIO::operator=(const PNGIO& source) {
130
	if (this != &source) {
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
		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
	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();
175
	}
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
	unsigned int w,h;
191
192
193
	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
	if (grid_w==0 || grid_h==0) {
214
215
216
217
		return 1.;
	}

	double min_factor = IOUtils::nodata;
218
	if (min_w!=IOUtils::unodata) { //min_w & min_w are read together
219
220
221
222
		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
	double max_factor = IOUtils::nodata;
225
	if (max_w!=IOUtils::unodata) { //max_w & max_h are read together
226
227
228
229
230
		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);
	}

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

	//only one size specification provided -> return its matching factor
237
	if (min_factor!=IOUtils::nodata)
238
239
240
		return min_factor;
	else
		return max_factor;
241
242
}

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

249
250
Grid2DObject PNGIO::scaleGrid(const Grid2DObject& grid_in)
{ //scale input image
251
	const double factor = getScaleFactor(grid_in.getNx(), grid_in.getNy());
252
	if (scaling=="nearest")
253
		return ResamplingAlgorithms2D::NearestNeighbour(grid_in, factor);
254
	else if (scaling=="bilinear")
255
		return ResamplingAlgorithms2D::BilinearResampling(grid_in, factor);
256
	else {
257
		ostringstream ss;
258
259
260
		ss << "Grid scaling algorithm \"" << scaling << "\" unknown";
		throw UnknownValueException(ss.str(), AT);
	}
261
}
262

263
void PNGIO::setFile(const std::string& filename, png_structp& png_ptr, png_infop& info_ptr, const size_t &width, const size_t &height)
264
{
265
	// Open file for writing (binary mode)
266
	if (!IOUtils::validFileAndPath(filename)) throw InvalidNameException(filename, AT);
267
	errno=0;
268
269
	fp = fopen(filename.c_str(), "wb");
	if (fp == NULL) {
270
		ostringstream ss;
271
		ss << "Error opening file \"" << filename << "\", possible reason: " << strerror(errno);
272
		throw AccessException(ss.str(), AT);
273
274
275
276
277
	}

	// Initialize write structure
	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
	if (png_ptr == NULL) {
278
		fclose(fp); fp=NULL;
279
		throw IOException("Could not allocate write structure. Are you running with the same version of libpng that you linked to?", AT);
280
281
282
283
284
	}

	// Initialize info structure
	info_ptr = png_create_info_struct(png_ptr);
	if (info_ptr == NULL) {
285
		fclose(fp); fp=NULL;
286
		png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
287
		free(png_ptr);
288
289
290
291
		throw IOException("Could not allocate info structure", AT);
	}

	// Setup Exception handling
292
#ifdef _MSC_VER
293
294
	#pragma warning(disable:4611) //the setjmp of libpng has been set up so that it can safely be called from c++
#endif
295
	if (setjmp(png_jmpbuf(png_ptr))) {
296
		closePNG(png_ptr, info_ptr, NULL);
297
		throw IOException("Error during png creation. Can not set jump pointer (I have no clue what it means too!)", AT);
298
299
300
301
	}

	png_init_io(png_ptr, fp);

302
	if (optimize_for_speed) png_set_compression_level(png_ptr, Z_BEST_SPEED);
303
304
	else png_set_compression_level(png_ptr, Z_BEST_COMPRESSION);

305
	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...
306
	if (indexed_png) png_set_compression_strategy(png_ptr, Z_RLE); //Z_DEFAULT_STRATEGY, Z_FILTERED, Z_HUFFMAN_ONLY, Z_RLE
307
308

	// Write header (8 bit colour depth). Full alpha channel with PNG_COLOR_TYPE_RGB_ALPHA
309
	if (indexed_png) {
310
		png_set_IHDR(png_ptr, info_ptr, static_cast<png_uint_32>(width), static_cast<png_uint_32>(height),
311
312
313
314
315
316
			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 {
317
		png_set_IHDR(png_ptr, info_ptr, static_cast<png_uint_32>(width), static_cast<png_uint_32>(height),
318
319
320
321
322
323
324
			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);
	}

325
326
327
	//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);
328
329
}

330
size_t PNGIO::setLegend(const size_t &ncols, const size_t &nrows, const double &min, const double &max, Array2D<double> &legend_array)
331
{
332
	if (has_legend) {
333
		const legend leg(static_cast<unsigned int>(nrows), min, max);
334
		legend_array = leg.getLegend();
335
		const size_t nx = legend_array.getNx();
336
337
338
		return (ncols+nx);
	} else {
		return ncols;
339
	}
340
341
}

342
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)
343
{
344
345
	const size_t ncols = grid.getNx();
	const size_t nrows = grid.getNy();
346

347
	// Allocate memory for one row (3 bytes per pixel - RGB)
348
	const unsigned char channels = (indexed_png)? 1 : 3; //4 for rgba
349
	png_bytep row = (png_bytep)calloc(channels*sizeof(png_byte), full_width);
350
	if (row==NULL) {
351
352
353
		fclose(fp); fp=NULL;
		png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
		free(png_ptr);
354
		throw IOException("Can not allocate row memory for PNG!", AT);
355
	}
356
357

	// Write image data
358
359
	if (indexed_png) {
		for (size_t y=nrows ; y-- > 0; ) {
360
			size_t x=0;
361
			for (; x<ncols ; x++) {
362
				const size_t i=x*channels;
363
				unsigned char index;
364
				gradient.getColor(grid(x,y), index);
365
				row[i]=static_cast<png_byte>(index);
366
			}
367
			for (; x<full_width; x++) {
368
				const size_t i=x*channels;
369
				unsigned char index;
370
				gradient.getColor(legend_array(x-ncols,y), index);
371
				row[i]=static_cast<png_byte>(index);
372
			}
373
			png_write_row(png_ptr, row);
374
		}
375
	} else {
376
		for (size_t y=nrows ; y -- > 0; ) {
377
			size_t x=0;
378
			for (; x<ncols ; x++) {
379
				const size_t i=x*channels;
380
381
382
				unsigned char r,g,b;
				bool a;
				gradient.getColor(grid(x,y), r,g,b,a);
383
				if (a==true) {
384
					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);
385
				} else {
386
					row[i]=static_cast<png_byte>(r); row[i+1]=static_cast<png_byte>(g); row[i+2]=static_cast<png_byte>(b);
387
388
				}
			}
389
			for (; x<full_width; x++) {
390
				const size_t i=x*channels;
391
392
393
				unsigned char r,g,b;
				bool a;
				gradient.getColor(legend_array(x-ncols,y), r,g,b,a);
394
				if (a==true) {
395
					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);
396
				} else {
397
					row[i]=static_cast<png_byte>(r); row[i+1]=static_cast<png_byte>(g); row[i+2]=static_cast<png_byte>(b);
398
				}
399
			}
400
			png_write_row(png_ptr, row);
401
402
		}
	}
403

404
405
	png_write_flush(png_ptr);
	png_free(png_ptr, row);
406
407
}

408
void PNGIO::setPalette(const Gradient &gradient, png_structp& png_ptr, png_infop& info_ptr, png_color *palette)
409
{
410
411
412
	std::vector<unsigned char> pal;
	size_t nr_colors;
	gradient.getPalette(pal, nr_colors);
413
	palette = (png_color*)calloc(sizeof (png_color), nr_colors); //ie: three png_bytes, each being an unsigned char
414
	for (size_t ii=0; ii<nr_colors; ii++) {
415
416
417
418
		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]);
419
	}
420
	png_set_PLTE(png_ptr, info_ptr, palette, static_cast<int>(nr_colors));
421
422
}

423
void PNGIO::closePNG(png_structp& png_ptr, png_infop& info_ptr, png_color *palette)
424
{
425
	png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
426
	if (indexed_png && palette!=NULL) free(palette);
427
	png_destroy_write_struct(&png_ptr, &info_ptr);
428
	fclose(fp); fp=NULL;
429
430
	free(info_ptr);
	free(png_ptr);
431
432
}

433
void PNGIO::write2DGrid(const Grid2DObject& grid_in, const std::string& filename)
434
{
435
	const string full_name = grid2dpath+"/"+filename;
436
	fp=NULL;
437
	png_color *palette=NULL;
438
439
440
441
442
	png_structp png_ptr=NULL;
	png_infop info_ptr=NULL;

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

446
447
448
449
	const double min = grid.grid2D.getMin();
	const double max = grid.grid2D.getMax();

	Gradient gradient(Gradient::heat, min, max, autoscale);
450
	if (indexed_png) gradient.setNrOfLevels(nr_levels);
451
452

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

455
	setFile(full_name, png_ptr, info_ptr, full_width, nrows);
456
457
	if (indexed_png) setPalette(gradient, png_ptr, info_ptr, palette);
	if (has_world_file) writeWorldFile(grid, full_name);
458
459
460
461
462
463

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

464
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
465
	png_write_end(png_ptr, NULL);
466

467
	closePNG(png_ptr, info_ptr, palette);
468
469
}

470
void PNGIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameters& parameter, const Date& date)
471
{
472
	std::string filename;
473
	if (parameter==MeteoGrids::DEM || parameter==MeteoGrids::SLOPE || parameter==MeteoGrids::AZI)
474
		filename = grid2dpath + "/" + MeteoGrids::getParameterName(parameter) + ".png";
475
476
477
478
479
	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";
	}
480

481
	fp=NULL;
482
	png_color *palette=NULL;
483
484
485
486
487
	png_structp png_ptr=NULL;
	png_infop info_ptr=NULL;

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

491
492
493
494
	double min = grid.grid2D.getMin();
	double max = grid.grid2D.getMax();

	Gradient gradient;
495
496
	if (parameter==MeteoGrids::DEM) {
		if (!autoscale) {
497
			min = 0.; //we want a 3000 snow line with a full scale legend
498
499
			max = 3500.;
			gradient.set(Gradient::terrain, min, max, autoscale); //max used as snow line reference
500
501
		} else
			gradient.set(Gradient::terrain, min, max, autoscale);
502
	} else if (parameter==MeteoGrids::SLOPE) {
503
		gradient.set(Gradient::slope, min, max, autoscale);
504
505
	} else if (parameter==MeteoGrids::SHADE) {
		if (!autoscale) {
506
507
508
			min = 0.; max = 1.;
		}
		gradient.set(Gradient::blktowhite, min, max, autoscale);
509
510
	} else if (parameter==MeteoGrids::AZI) {
		if (!autoscale) {
511
512
513
			min = 0.;
			max = 360.;
		}
514
		gradient.set(Gradient::azi, min, max, autoscale);
515
516
	} else if (parameter==MeteoGrids::DW) {
		if (!autoscale) {
517
518
519
520
			min = 0.;
			max = 360.;
		}
		gradient.set(Gradient::azi, min, max, autoscale);
521
522
	} else if (parameter==MeteoGrids::HS) {
		if (!autoscale) {
523
			min = 0.; max = 2.5;
524
525
		}
		gradient.set(Gradient::blue, min, max, autoscale);
526
	} else if (parameter==MeteoGrids::TA) {
527
		grid.grid2D -= Cst::t_water_freezing_pt; //convert to celsius
528
		if (!autoscale) {
529
			min = -15.; max = 15.;
530
531
532
533
		} else {
			min -= Cst::t_water_freezing_pt;
			max -= Cst::t_water_freezing_pt;
		}
534
		gradient.set(Gradient::heat, min, max, autoscale);
535
	} else if (parameter==MeteoGrids::TSS) {
536
		grid.grid2D -= Cst::t_water_freezing_pt; //convert to celsius
537
		if (!autoscale) {
538
			min = -20.; max = 5.;
539
540
541
542
543
		} else {
			min -= Cst::t_water_freezing_pt;
			max -= Cst::t_water_freezing_pt;
		}
		gradient.set(Gradient::freeze, min, max, autoscale);
544
545
546
547
548
549
550
551
552
	} else if (parameter==MeteoGrids::TSOIL) {
		grid.grid2D -= Cst::t_water_freezing_pt; //convert to celsius
		if (!autoscale) {
			min = -5.; max = 5.;
		} else {
			min -= Cst::t_water_freezing_pt;
			max -= Cst::t_water_freezing_pt;
		}
		gradient.set(Gradient::heat, min, max, autoscale);
553
554
	} else if (parameter==MeteoGrids::RH) {
		if (!autoscale) {
555
556
			min = 0.; max = 1.;
		}
557
		gradient.set(Gradient::bg_isomorphic, min, max, autoscale);
558
559
	} else if (parameter==MeteoGrids::P) {
		if (!autoscale) {
560
561
562
563
564
565
566
567
568
			//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);
		}
569
570
	} else if (parameter==MeteoGrids::ALB) {
		if (!autoscale) {
571
572
			min = 0.; max = 1.;
		}
573
		gradient.set(Gradient::blktowhite, min, max, autoscale);
574
575
	} else if (parameter==MeteoGrids::TAU_CLD) {
		if (!autoscale) {
576
577
578
			min = 0.; max = 1.;
		}
		gradient.set(Gradient::blktowhite, min, max, autoscale);
579
580
	} else if (parameter==MeteoGrids::ISWR) {
		if (!autoscale) {
581
			min = 0.; max = 800.;
582
583
		}
		gradient.set(Gradient::heat, min, max, autoscale);
584
585
	} else if (parameter==MeteoGrids::ILWR) {
		if (!autoscale) {
586
587
588
			min = 200.; max = 500.;
		}
		gradient.set(Gradient::heat, min, max, autoscale);
589
590
	} else if (parameter==MeteoGrids::SWE) {
		if (!autoscale) {
591
			min = 0.; max = 250.;
592
		}
593
		gradient.set(Gradient::blue_pink, min, max, autoscale);
594
595
596
	} else {
		gradient.set(Gradient::heat, min, max, autoscale);
	}
597
	gradient.setNrOfLevels(nr_levels);
598
599

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

602
	setFile(filename, png_ptr, info_ptr, full_width, nrows);
603
604
	if (indexed_png) setPalette(gradient, png_ptr, info_ptr, palette);
	if (has_world_file) writeWorldFile(grid, filename);
605
606
607
608

	createMetadata(grid);
	metadata_key.push_back("Title"); //adding title
	metadata_text.push_back( MeteoGrids::getParameterName(parameter)+" on "+date.toString(Date::ISO) );
609
610
	metadata_key.push_back("Simulation Date");
	metadata_text.push_back( date.toString(Date::ISO) );
611
612
	metadata_key.push_back("Simulation Parameter");
	metadata_text.push_back( MeteoGrids::getParameterName(parameter) );
613
614
	writeMetadata(png_ptr, info_ptr);

615
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
616
	png_write_end(png_ptr, NULL);
617

618
	closePNG(png_ptr, info_ptr, palette);
619
620
621
622
623
624
}

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;
625
	Coords world_ref = grid_in.llcorner;
626
	world_ref.setProj(coordout, coordoutparam);
627
	world_ref.moveByXY(.5*cellsize, (double(grid_in.getNy())+.5)*cellsize); //moving to center of upper left cell
628

629
	if (!IOUtils::validFileAndPath(world_file)) throw InvalidNameException(world_file, AT);
630
	std::ofstream fout(world_file.c_str(), ios::out);
631
	if (fout.fail()) {
632
		throw AccessException(world_file, AT);
633
634
635
636
637
638
639
640
641
642
643
	}

	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();
644
		throw AccessException("Failed when writing to PNG world file \""+world_file+"\"", AT);
645
646
647
	}

	fout.close();
648
649
}

650
void PNGIO::createMetadata(const Grid2DObject& grid)
651
{
652
653
	const double lat = grid.llcorner.getLat();
	const double lon = grid.llcorner.getLon();
654
	ostringstream ss;
655

656
657
658
	metadata_key.clear();
	metadata_text.clear();

659
660
661
662
	metadata_key.push_back("Creation Time");
	Date cr_date;
	cr_date.setFromSys();
	metadata_text.push_back( cr_date.toString(Date::ISO) );
663
664
665
666
667
668
669
670
671
	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());
672
673
	metadata_key.push_back("LL_Latitude");
	ss.str(""); ss << fixed << setprecision(7) << lat;
674
	metadata_text.push_back(ss.str());
675
676
677
678
679
680
681
682
683
684
685
	metadata_key.push_back("LL_Longitude");
	ss.str(""); ss << fixed << setprecision(7) << lon;
	metadata_text.push_back(ss.str());
	
	Coords UR(grid.llcorner);
	UR.moveByXY( grid.cellsize*static_cast<double>(grid.getNx()) , grid.cellsize*static_cast<double>(grid.getNy()) );
	metadata_key.push_back("UR_Latitude");
	ss.str(""); ss << fixed << setprecision(7) << UR.getLat();
	metadata_text.push_back(ss.str());
	metadata_key.push_back("UR_Longitude");
	ss.str(""); ss << fixed << setprecision(7) << UR.getLon();
686
	metadata_text.push_back(ss.str());
687

688
	if (lat<0.) {
689
690
691
692
693
694
695
696
697
698
		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));
	}
699
	if (lon<0.) {
700
701
702
703
704
705
706
707
708
709
710
711
712
713
		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)
{
714
	const size_t max_len = 79; //according to the official specs' recommendation
715
	const size_t nr = metadata_key.size();
716
717
718
	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);
719

720
	for (size_t ii=0; ii<nr; ii++) {
721
722
		key[ii] = (char *)calloc(sizeof(char), max_len);
		text[ii] = (char *)calloc(sizeof(char), max_len);
723
		strncpy(key[ii], metadata_key[ii].c_str(), max_len-1); //in case the '\0' was not counted by maxlen
724
		key[ii][max_len-1] = '\0'; //make sure it is null terminated
725
		strncpy(text[ii], metadata_text[ii].c_str(), max_len-1);
726
727
728
		text[ii][max_len-1] = '\0'; //make sure it is null terminated
		info_text[ii].key = key[ii];
		info_text[ii].text = text[ii];
729
730
		info_text[ii].compression = PNG_TEXT_COMPRESSION_NONE;
	}
731

732
	png_set_text(png_ptr, info_ptr, info_text, static_cast<int>(nr));
733
734
	png_write_info(png_ptr, info_ptr);

735
	free(info_text);
736
	for (size_t ii=0; ii<nr; ii++) {
737
738
739
740
741
742
743
744
745
746
747
748
		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;

749
	std::ostringstream dms;
750
751
	dms << d << "/1 " << static_cast<int>(m*100) << "/100 " << fixed << setprecision(6) << s << "/1";
	return dms.str();
752
753
754
}

} //namespace