WSL/SLF GitLab Repository

PNGIO.cc 28.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  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/>.
*/
18
#include <meteoio/plugins/PNGIO.h>
19
#include <meteoio/ResamplingAlgorithms2D.h>
20
#include <meteoio/Graphics.h>
21
#include <meteoio/FileUtils.h>
22
#include <meteoio/meteoLaws/Meteoconst.h>
23

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

using namespace std;

namespace mio {
/**
34
 * @page pngio PNGIO
35
 * @section pngio_format Format
36
37
38
39
 * 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).
40
 * If a grid containing no data (ie: size 0x0) is sent to the plugin, then no file will be written.
41
 * Finally, the naming scheme for meteo grids should be: YYYY-MM-DDTHH.mm_{MeteoGrids::Parameters}.png
42
 *
43
 * @section pngio_units Units
44
 * All units are MKSA except temperatures that are expressed in celcius.
45
 *
46
 * @section pngio_keywords Keywords
47
 * This plugin uses the following keywords:
48
49
 * - COORDSYS: input coordinate system (see Coords) specified in the [Output] section
 * - COORDPARAM: extra input coordinates parameters (see Coords) specified in the [Output] section
50
51
52
53
54
55
56
 * - 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)
57
58
 *
 * Advanced parameters (ie: don't mess up with them if you don't know what you're doing):
59
 * - PNG_INDEXED: create an indexed PNG? (default=true)
60
 * - PNG_NR_LEVELS: number of colors to use (less=smaller files, but it must be at least 5 and less than 255. default=30)
61
 * - PNG_SPEED_OPTIMIZE: optimize file creation for speed? (default=true, otherwise optimize for file size)
62
63
 *
 * 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.
64
65
66
67
68
69
70
71
72
73
74
 * 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
75
 *
76
 * @section pngio_example Example use
77
78
79
80
81
82
 * @code
 * GRID2D = PNG
 * png_legend = false
 * png_min_size = 400x400
 * png_max_size = 1366*768
 * @endcode
83
 *
84
 * @section pngio_compilation
85
86
 * 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.
 *
87
88
89
90
91
 * 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>).
 * 
92
93
94
95
96
97
98
99
100
101
 * 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
102
103
104
 */

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

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

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

131
PNGIO& PNGIO::operator=(const PNGIO& source) {
132
	if (this != &source) {
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
		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;
}

Mathias Bavay's avatar
Mathias Bavay committed
154
155
156
157
158
159
PNGIO::~PNGIO() throw() 
{
	if (fp!=NULL) fclose(fp); 
	fp=NULL;
}

160
161
void PNGIO::setOptions()
{
162
	//default values have been set by the constructors
163
	cfg.getValue("COORDSYS", "Output", coordout);
164
	cfg.getValue("COORDPARAM", "Output", coordoutparam, IOUtils::nothrow);
165
	cfg.getValue("GRID2DPATH", "Output", grid2dpath);
166
	//cfg.getValue("TIME_ZONE", "Output", tz_out, IOUtils::nothrow);
167

168
169
	//get size specifications
	std::string min_size, max_size;
170
	cfg.getValue("PNG_MIN_SIZE", "Output", min_size, IOUtils::nothrow);
171
	if (!min_size.empty()) parse_size(min_size, min_w, min_h);
172
	cfg.getValue("PNG_MAX_SIZE", "Output", max_size, IOUtils::nothrow);
173
	if (!max_size.empty()) parse_size(max_size, max_w, max_h);
174

175
176
177
178
	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);
179

180
	if (has_legend) { //we need to save room for the legend
181
182
		if (min_w!=IOUtils::unodata) min_w -= Legend::getLegendWidth();
		if (max_w!=IOUtils::unodata) max_w -= Legend::getLegendWidth();
183
	}
184

185
186
	cfg.getValue("PNG_INDEXED", "Output", indexed_png, IOUtils::nothrow);
	cfg.getValue("PNG_SPEED_OPTIMIZE", "Output", optimize_for_speed, IOUtils::nothrow);
187
	unsigned int tmp=IOUtils::unodata;
188
	cfg.getValue("PNG_NR_LEVELS", "Output", tmp, IOUtils::nothrow);
189
	if (tmp!=IOUtils::unodata && (tmp>255 || tmp<5)) {
190
191
		throw InvalidFormatException("PNG_NR_LEVELS must be between 5 and 255!", AT);
	}
192
	if (tmp!=IOUtils::unodata) nr_levels=static_cast<unsigned char>(tmp);
193
194
}

195
void PNGIO::parse_size(const std::string& size_spec, size_t& width, size_t& height)
196
197
{
	char rest[32] = "";
198
	unsigned int w,h;
199
200
201
	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) {
202
		std::ostringstream ss;
203
204
205
		ss << "Can not parse PNGIO size specification \"" << size_spec << "\"";
		throw InvalidFormatException(ss.str(), AT);
	}
206

Mathias Bavay's avatar
Mathias Bavay committed
207
208
	width = static_cast<size_t>( w );
	height = static_cast<size_t>( h );
209

210
211
212
	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
Mathias Bavay's avatar
Mathias Bavay committed
213
		throw InvalidFormatException("Invalid PNGIO size specification \"" + size_spec + "\"", AT);
214
215
216
	}
}

Mathias Bavay's avatar
Mathias Bavay committed
217
double PNGIO::getScaleFactor(const size_t& grid_w, const size_t& grid_h) const
218
{
219
	if (grid_w==0 || grid_h==0) {
220
221
222
223
		return 1.;
	}

	double min_factor = IOUtils::nodata;
224
	if (min_w!=IOUtils::unodata) { //min_w & min_w are read together
225
226
227
228
		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);
	}
229

230
	double max_factor = IOUtils::nodata;
231
	if (max_w!=IOUtils::unodata) { //max_w & max_h are read together
232
233
234
235
236
		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);
	}

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

	//only one size specification provided -> return its matching factor
243
	if (min_factor!=IOUtils::nodata)
244
245
246
		return min_factor;
	else
		return max_factor;
247
248
}

Mathias Bavay's avatar
Mathias Bavay committed
249
Grid2DObject PNGIO::scaleGrid(const Grid2DObject& grid_in) const
250
{ //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 (!FileUtils::validFileAndPath(filename)) throw InvalidNameException(filename, AT);
267
	errno=0;
268
	fp = fopen(filename.c_str(), "wb");
269
270
	if (fp == NULL)
		throw AccessException("Error opening file \""+filename+"\", possible reason: "+std::string(strerror(errno)), AT);
271
272
273
274

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

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

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

	png_init_io(png_ptr, fp);

299
	if (optimize_for_speed) png_set_compression_level(png_ptr, Z_BEST_SPEED);
300
301
	else png_set_compression_level(png_ptr, Z_BEST_COMPRESSION);

302
	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...
303
	if (indexed_png) png_set_compression_strategy(png_ptr, Z_RLE); //Z_DEFAULT_STRATEGY, Z_FILTERED, Z_HUFFMAN_ONLY, Z_RLE
304
305

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

322
323
324
	//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);
325
326
}

Mathias Bavay's avatar
Mathias Bavay committed
327
size_t PNGIO::setLegend(const size_t &ncols, const size_t &nrows, const double &min, const double &max, Array2D<double> &legend_array) const
328
{
329
	if (has_legend) {
330
331
		const Legend legend(static_cast<unsigned int>(nrows), min, max);
		legend_array = legend.getLegend();
332
		const size_t nx = legend_array.getNx();
333
334
335
		return (ncols+nx);
	} else {
		return ncols;
336
	}
337
338
}

Mathias Bavay's avatar
Mathias Bavay committed
339
340
341
342
343
void PNGIO::setPalette(const Gradient &gradient, png_structp& png_ptr, png_infop& info_ptr, png_color *palette)
{
	std::vector<unsigned char> pal;
	size_t nr_colors;
	gradient.getPalette(pal, nr_colors);
344
	palette = (png_color*)calloc(nr_colors, sizeof (png_color)); //ie: three png_bytes, each being an unsigned char
Mathias Bavay's avatar
Mathias Bavay committed
345
346
347
348
349
350
351
352
353
	for (size_t ii=0; ii<nr_colors; ii++) {
		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]);
	}
	png_set_PLTE(png_ptr, info_ptr, palette, static_cast<int>(nr_colors));
}

354
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)
355
{
356
357
	const size_t ncols = grid.getNx();
	const size_t nrows = grid.getNy();
358

359
	// Allocate memory for one row (3 bytes per pixel - RGB)
360
	const unsigned char channels = (indexed_png)? 1 : 3; //4 for rgba
361
	png_bytep row = (png_bytep)calloc(full_width, channels*sizeof(png_byte));
362
	if (row==NULL) {
363
364
365
		fclose(fp); fp=NULL;
		png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
		free(png_ptr);
366
		throw IOException("Can not allocate row memory for PNG!", AT);
367
	}
368
369

	// Write image data
370
371
	if (indexed_png) {
		for (size_t y=nrows ; y-- > 0; ) {
372
			size_t x=0;
373
			for (; x<ncols ; x++) {
374
				const size_t i=x*channels;
375
				unsigned char index;
376
				gradient.getColor(grid(x,y), index);
377
				row[i]=static_cast<png_byte>(index);
378
			}
379
			for (; x<full_width; x++) {
380
				const size_t i=x*channels;
381
				unsigned char index;
382
				gradient.getColor(legend_array(x-ncols,y), index);
383
				row[i]=static_cast<png_byte>(index);
384
			}
385
			png_write_row(png_ptr, row);
386
		}
387
	} else {
388
		for (size_t y=nrows ; y -- > 0; ) {
389
			size_t x=0;
390
			for (; x<ncols ; x++) {
391
				const size_t i=x*channels;
392
393
394
				unsigned char r,g,b;
				bool a;
				gradient.getColor(grid(x,y), r,g,b,a);
395
				if (a==true) {
396
					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);
397
				} else {
398
					row[i]=static_cast<png_byte>(r); row[i+1]=static_cast<png_byte>(g); row[i+2]=static_cast<png_byte>(b);
399
400
				}
			}
401
			for (; x<full_width; x++) {
402
				const size_t i=x*channels;
403
404
405
				unsigned char r,g,b;
				bool a;
				gradient.getColor(legend_array(x-ncols,y), r,g,b,a);
406
				if (a==true) {
407
					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);
408
				} else {
409
					row[i]=static_cast<png_byte>(r); row[i+1]=static_cast<png_byte>(g); row[i+2]=static_cast<png_byte>(b);
410
				}
411
			}
412
			png_write_row(png_ptr, row);
413
414
		}
	}
415

416
417
	png_write_flush(png_ptr);
	png_free(png_ptr, row);
418
419
}

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

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

	//scale input image
439
	const Grid2DObject grid( scaleGrid(grid_in) );
440
	const size_t ncols = grid.getNx(), nrows = grid.getNy();
441
	if (ncols==0 || nrows==0) return;
442

443
444
445
446
	const double min = grid.grid2D.getMin();
	const double max = grid.grid2D.getMax();

	Gradient gradient(Gradient::heat, min, max, autoscale);
447
	if (indexed_png) gradient.setNrOfLevels(nr_levels);
448
449

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

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

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

461
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
462
	png_write_end(png_ptr, NULL);
463

464
	closePNG(png_ptr, info_ptr, palette);
465
466
}

467
void PNGIO::write2DGrid(const Grid2DObject& grid_in, const MeteoGrids::Parameters& parameter, const Date& date)
468
{
469
470
471
	const bool isNormalParam = (parameter!=MeteoGrids::DEM && parameter!=MeteoGrids::SLOPE && parameter!=MeteoGrids::AZI);
	std::string date_str = (isNormalParam)? date.toString(Date::ISO)+"_" : "";
	std::replace( date_str.begin(), date_str.end(), ':', '.');
472
	const std::string filename( grid2dpath + "/" + date_str + MeteoGrids::getParameterName(parameter) + ".png" );
473

474
	fp=NULL;
475
	png_color *palette=NULL;
476
477
478
479
	png_structp png_ptr=NULL;
	png_infop info_ptr=NULL;

	//scale input image
480
	Grid2DObject grid( scaleGrid(grid_in) );
481
	const size_t ncols = grid.getNx(), nrows = grid.getNy();
482
	if (ncols==0 || nrows==0) return;
483

484
485
486
487
	double min = grid.grid2D.getMin();
	double max = grid.grid2D.getMax();

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

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

598
	setFile(filename, png_ptr, info_ptr, full_width, nrows);
599
600
	if (indexed_png) setPalette(gradient, png_ptr, info_ptr, palette);
	if (has_world_file) writeWorldFile(grid, filename);
601
602
603
604

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

611
	writeDataSection(grid, legend_array, gradient, full_width, png_ptr, info_ptr);
612
	png_write_end(png_ptr, NULL);
613

614
	closePNG(png_ptr, info_ptr, palette);
615
616
}

Mathias Bavay's avatar
Mathias Bavay committed
617
void PNGIO::writeWorldFile(const Grid2DObject& grid_in, const std::string& filename) const
618
{
619
	const std::string world_file( FileUtils::removeExtension(filename)+".pnw" );
620
	if (!FileUtils::validFileAndPath(world_file)) throw InvalidNameException(world_file, AT);
621
	std::ofstream fout(world_file.c_str(), ios::out);
Mathias Bavay's avatar
Mathias Bavay committed
622
	if (fout.fail()) throw AccessException(world_file, AT);
623

624
625
626
627
628
	const double cellsize = grid_in.cellsize;
	Coords world_ref( grid_in.llcorner );
	world_ref.setProj(coordout, coordoutparam);
	world_ref.moveByXY(.5*cellsize, (double(grid_in.getNy())+.5)*cellsize); //moving to center of upper left cell

629
630
631
632
633
634
635
636
637
	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();
638
		throw AccessException("Failed when writing to PNG world file \""+world_file+"\"", AT);
639
640
641
	}

	fout.close();
642
643
}

644
void PNGIO::createMetadata(const Grid2DObject& grid)
645
{
646
647
	const double lat = grid.llcorner.getLat();
	const double lon = grid.llcorner.getLon();
648
	ostringstream ss;
649

650
651
652
	metadata_key.clear();
	metadata_text.clear();

653
654
655
656
	metadata_key.push_back("Creation Time");
	Date cr_date;
	cr_date.setFromSys();
	metadata_text.push_back( cr_date.toString(Date::ISO) );
657
658
659
660
661
662
663
664
665
	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());
666
667
	metadata_key.push_back("LL_Latitude");
	ss.str(""); ss << fixed << setprecision(7) << lat;
668
	metadata_text.push_back(ss.str());
669
670
671
672
673
674
675
676
677
678
679
	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();
680
	metadata_text.push_back(ss.str());
681

682
	if (lat<0.) {
683
684
685
686
687
688
689
690
691
692
		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));
	}
693
	if (lon<0.) {
694
695
696
697
698
699
700
701
702
703
704
705
706
707
		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)
{
708
	const size_t max_len = 79; //according to the official specs' recommendation
709
	const size_t nr = metadata_key.size();
710
711
712
	png_text *info_text = (png_text *)calloc(nr, sizeof(png_text));
	char **key = (char**)calloc(nr, sizeof(char)*max_len);
	char **text = (char**)calloc(nr, sizeof(char)*max_len);
713

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

726
	png_set_text(png_ptr, info_ptr, info_text, static_cast<int>(nr));
727
728
	png_write_info(png_ptr, info_ptr);

729
	free(info_text);
730
	for (size_t ii=0; ii<nr; ii++) {
731
732
733
734
735
736
737
738
739
740
741
742
		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;

743
	std::ostringstream dms;
744
	dms << d << "/1 " << static_cast<int>(m*100.) << "/100 " << fixed << setprecision(6) << s << "/1";
745
	return dms.str();
746
747
748
}

} //namespace