WSL/SLF GitLab Repository

PNGIO.cc 29.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
	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
	// Open file for writing (binary mode)
321
	if (!IOUtils::validFileAndPath(filename)) throw InvalidFileNameException(filename, AT);
322
	errno=0;
323
324
	fp = fopen(filename.c_str(), "wb");
	if (fp == NULL) {
325
		ostringstream ss;
326
		ss << "Error opening file \"" << filename << "\", possible reason: " << strerror(errno);
327
		throw FileAccessException(ss.str(), AT);
328
329
330
331
332
	}

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

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

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

	png_init_io(png_ptr, fp);

357
	if (optimize_for_speed) png_set_compression_level(png_ptr, Z_BEST_SPEED);
358
359
	else png_set_compression_level(png_ptr, Z_BEST_COMPRESSION);

360
	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...
361
	if (indexed_png) png_set_compression_strategy(png_ptr, Z_RLE); //Z_DEFAULT_STRATEGY, Z_FILTERED, Z_HUFFMAN_ONLY, Z_RLE
362
363

	// Write header (8 bit colour depth). Full alpha channel with PNG_COLOR_TYPE_RGB_ALPHA
364
	if (indexed_png) {
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
		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);
	}

380
381
382
	//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);
383
384
}

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

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

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

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

459
460
	png_write_flush(png_ptr);
	png_free(png_ptr, row);
461
462
}

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

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

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

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

501
502
503
504
	const double min = grid.grid2D.getMin();
	const double max = grid.grid2D.getMax();

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

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

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

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

519
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
520
	png_write_end(png_ptr, NULL);
521

522
	closePNG(png_ptr, info_ptr, palette);
523
524
}

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

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

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

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

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

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

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

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

661
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
662
	png_write_end(png_ptr, NULL);
663

664
	closePNG(png_ptr, info_ptr, palette);
665
666
667
668
669
670
}

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

675
676
	if (!IOUtils::validFileAndPath(world_file)) throw InvalidFileNameException(world_file, AT);
	std::ofstream fout(world_file.c_str(), ios::out);
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
	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();
694
695
}

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

702
703
704
	metadata_key.clear();
	metadata_text.clear();

705
706
707
708
	metadata_key.push_back("Creation Time");
	Date cr_date;
	cr_date.setFromSys();
	metadata_text.push_back( cr_date.toString(Date::ISO) );
709
710
711
712
713
714
715
716
717
	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());
718
719
	metadata_key.push_back("LL_Latitude");
	ss.str(""); ss << fixed << setprecision(7) << lat;
720
	metadata_text.push_back(ss.str());
721
722
723
724
725
726
727
728
729
730
731
	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();
732
	metadata_text.push_back(ss.str());
733

734
	if (lat<0.) {
735
736
737
738
739
740
741
742
743
744
		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));
	}
745
	if (lon<0.) {
746
747
748
749
750
751
752
753
754
755
756
757
758
759
		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)
{
760
	const size_t max_len = 79; //according to the official specs' recommendation
761
	const size_t nr = metadata_key.size();
762
763
764
	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);
765

766
	for (size_t ii=0; ii<nr; ii++) {
767
768
		key[ii] = (char *)calloc(sizeof(char), max_len);
		text[ii] = (char *)calloc(sizeof(char), max_len);
769
770
771
772
		strncpy(key[ii], metadata_key[ii].c_str(), max_len-1); //in case the '\0' was not counted by maxlen
		strncpy(text[ii], metadata_text[ii].c_str(), max_len-1);
		info_text[ii].key = key[ii]+'\0'; //strncpy does not always '\0' terminate
		info_text[ii].text = text[ii]+'\0';
773
774
		info_text[ii].compression = PNG_TEXT_COMPRESSION_NONE;
	}
775

776
	png_set_text(png_ptr, info_ptr, info_text, static_cast<int>(nr));
777
778
	png_write_info(png_ptr, info_ptr);

779
	free(info_text);
780
	for (size_t ii=0; ii<nr; ii++) {
781
782
783
784
785
786
787
788
789
790
791
792
		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;

793
	std::ostringstream dms;
794
795
	dms << d << "/1 " << static_cast<int>(m*100) << "/100 " << fixed << setprecision(6) << s << "/1";
	return dms.str();
796
797
798
}

} //namespace