/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. 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 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_NETCDF_IO_H
#define _GCTL_NETCDF_IO_H

// get library configuration
#include "../gctl_config.h"
#ifdef GCTL_NETCDF

// netcdf cpp head file
#include "netcdfcxx_legacy/netcdfcpp.h"
#include "typeinfo"

#include "../core/array.h"
#include "../core/matrix.h"
#include "../utility/stream.h"

namespace gctl
{
	/**
	 * @brief      显示 Netcdf 格式文件信息
	 *
	 * @param[in]  filename    读入文件名（无后缀）
	 * @param      out_stream  输出流 默认为标准输出
	 */
	void show_netcdf_info(std::string filename, std::ostream& out_stream = std::cout);


	template <typename T>
	void read_netcdf_axis(std::string filename, array<T> &out_x, std::string xname)
	{
		if (filename == "")
		{
			throw domain_error("The input filename is empty. From read_netcdf_axis(...)");
		}

		// 判断文件名结尾是否为.nc如果不是则添加.nc结尾
		std::string full_name;
		if (filename.length() <= 3 || filename.substr(filename.length()-3, filename.length()) != ".nc")
		{
			full_name = filename + ".nc";
		}
		else full_name = filename;

		NcFile input_nc(full_name.c_str(), NcFile::ReadOnly);

		if (!input_nc.is_valid())
		{
			throw runtime_error("Can't open file: " + full_name + ". From read_netcdf_axis(...)");
		}

		NcVar *x_var = input_nc.get_var(xname.c_str());
		if (!x_var->is_valid())
		{
			throw runtime_error("Variable not found: " + xname + ". From read_netcdf_axis(...)");
		}

		NcType axis_type = x_var->type();
		if ((axis_type == ncByte   && typeid(T).name() != typeid(signed char).name()) || 
			(axis_type == ncChar   && typeid(T).name() != typeid(char).name()) || 
			(axis_type == ncShort  && typeid(T).name() != typeid(short).name()) || 
			(axis_type == ncInt    && typeid(T).name() != typeid(int).name()) || 
			(axis_type == ncFloat  && typeid(T).name() != typeid(float).name()) || 
			(axis_type == ncDouble && typeid(T).name() != typeid(double).name()))
		{
			throw runtime_error("Incompatible extracting data type for the axis: " + xname + ". From read_netcdf_axis(...)");
		}

		if (!out_x.empty()) out_x.clear();
		out_x.resize(x_var->num_vals());
		x_var->get(out_x.get(), out_x.size());
		return;
	}

	/**
	 * @brief      读入 Netcdf 格式的规则网格数据文件
	 * 
	 * @warning    目前只支持 Classic 格式的 .nc 文件。
	 *
	 * @param      out_data  返回数组（一维数据，左下到右上，行优先排序）的指针，默认为 nullptr
	 * @param[in]  filename  读入文件名（无后缀）
	 * @param[in]  xname     文件中 x 坐标轴的变量名，默认为 x
	 * @param[in]  yname     文件中 y 坐标轴的变量名，默认为 y
	 * @param[in]  dataname  文件中网格数据的变量名，默认为 z
	 */
	template <typename T>
	void read_netcdf_grid(std::string filename, array<T> &out_data, std::string xname = "x", std::string yname = "y", std::string dataname = "z")
	{
		if (filename == "")
		{
			throw domain_error("The input filename is empty. From read_netcdf_grid(...)");
		}

		// 判断文件名结尾是否为.nc如果不是则添加.nc结尾
		std::string full_name;
		if (filename.length() <= 3 || filename.substr(filename.length()-3, filename.length()) != ".nc")
		{
			full_name = filename + ".nc";
		}
		else full_name = filename;

		NcFile input_nc(full_name.c_str(), NcFile::ReadOnly);

		if (!input_nc.is_valid())
		{
			throw runtime_error("Can't open file: " + filename + ".nc . From read_netcdf_grid(...)");
		}

		NcVar *x_var = input_nc.get_var(xname.c_str());
		if (!x_var->is_valid())
		{
			throw runtime_error("Variable not found: " + xname + ". From read_netcdf_grid(...)");
		}

		NcVar *y_var = input_nc.get_var(yname.c_str());
		if (!y_var->is_valid())
		{
			throw runtime_error("Variable not found: " + yname + ". From read_netcdf_grid(...)");
		}

		NcVar *z_var = input_nc.get_var(dataname.c_str());
		if (!z_var->is_valid())
		{
			throw runtime_error("Variable not found: " + dataname + ". From read_netcdf_grid(...)");
		}
/*
		NcType val_type = z_var->type();
		if ((val_type == ncByte   && typeid(T).name() != typeid(signed char).name()) || 
			(val_type == ncChar   && typeid(T).name() != typeid(char).name()) || 
			(val_type == ncShort  && typeid(T).name() != typeid(short).name()) || 
			(val_type == ncInt    && typeid(T).name() != typeid(int).name()) || 
			(val_type == ncFloat  && typeid(T).name() != typeid(float).name()) || 
			(val_type == ncDouble && typeid(T).name() != typeid(double).name()))
		{
			throw runtime_error("Incompatible extracting data type for the data: " + dataname + ". From read_netcdf_grid(...)");
		}
*/
		if (!out_data.empty()) out_data.clear();
		out_data.resize(z_var->num_vals());

		if (z_var->type() == ncFloat)
		{
			array<float> ex_data(z_var->num_vals());
			z_var->get(ex_data.get(), y_var->num_vals(), x_var->num_vals());

			for (size_t i = 0; i < out_data.size(); i++)
			{
				out_data[i] = (T) ex_data[i];
			}
		}
		else if (z_var->type() == ncDouble)
		{
			array<double> ex_data(z_var->num_vals());
			z_var->get(ex_data.get(), y_var->num_vals(), x_var->num_vals());

			for (size_t i = 0; i < out_data.size(); i++)
			{
				out_data[i] = (T) ex_data[i];
			}
		}
		else
		{
			throw runtime_error("Unsupported data type(s) yet. From read_netcdf_grid(...)");
		}
		return;
	}

	/**
	 * @brief      读入 Netcdf 格式的规则网格的所有数据文件
	 * 
	 * @warning    目前只支持 Classic 格式的 .nc 文件。
	 *
	 * @param      out_data  返回数组（一维数据，左下到右上，行优先排序）的指针，默认为 nullptr
	 * @param[in]  filename  读入文件名（无后缀）
	 * @param[in]  xname     文件中 x 坐标轴的变量名，默认为 x
	 * @param[in]  yname     文件中 y 坐标轴的变量名，默认为 y
	 */
	template <typename T>
	void read_netcdf_grid(std::string filename, matrix<T> &out_data, array<std::string> &out_name, std::string xname = "x", std::string yname = "y")
	{
		if (filename == "")
		{
			throw domain_error("The input filename is empty. From read_netcdf_grid(...)");
		}

		// 判断文件名结尾是否为.nc如果不是则添加.nc结尾
		std::string full_name;
		if (filename.length() <= 3 || filename.substr(filename.find_last_of('.')) != ".nc")
		{
			full_name = filename + ".nc";
		}
		else full_name = filename;

		NcFile input_nc(full_name.c_str(), NcFile::ReadOnly);

		if (!input_nc.is_valid())
		{
			throw runtime_error("Can't open file: " + filename + ".nc . From read_netcdf_grid(...)");
		}

		int row = 0, col = 0;
		NcVar *a_var;
		NcDim *xdim, *ydim;
		NcValues *var_val;
		for (int n = 0; n < input_nc.num_vars(); n++)
		{
			a_var = input_nc.get_var(n);

			if (a_var->num_dims() == 2)
			{
				xdim = a_var->get_dim(1);
				ydim = a_var->get_dim(0);
				if (xdim->name() == xname && ydim->name() == yname)
				{
					row++;
					col = a_var->num_vals();
				}
			}
		}

		out_name.resize(row);
		out_data.resize(row, col);

		int c = 0;
		for (int n = 0; n < input_nc.num_vars(); n++)
		{
			a_var = input_nc.get_var(n);

			if (a_var->num_dims() == 2)
			{
				xdim = a_var->get_dim(1);
				ydim = a_var->get_dim(0);
				if (xdim->name() == xname && ydim->name() == yname)
				{
					out_name[c] = a_var->name();
					var_val = a_var->values();
					for (int i = 0; i < a_var->num_vals(); i++)
					{
						if (a_var->type() == ncFloat) out_data[c][i] = var_val->as_float(i);
						if (a_var->type() == ncDouble) out_data[c][i] = var_val->as_double(i);
					}
					c++;
				}
			}
		}
		return;
	}

	/**
	 * @brief      读入 Netcdf 格式的规则网格数据文件
	 * 
	 * @warning    目前只支持 Classic 格式的 .nc 文件。
	 *
	 * @param      out_x     返回一维数组
	 * @param      out_y     返回一维数组
	 * @param      out_data  返回二维数组
	 * @param[in]  filename  读入文件名（无后缀）
	 * @param[in]  xname     文件中 x 坐标轴的变量名，默认为 x
	 * @param[in]  yname     文件中 y 坐标轴的变量名，默认为 y
	 * @param[in]  dataname  文件中网格数据的变量名，默认为 z
	 */
	template <typename T, typename P>
	void read_netcdf_grid(std::string filename, array<P> &out_x, array<P> &out_y, matrix<T> &out_data, std::string xname = "x", std::string yname = "y", std::string dataname = "z")
	{
		array<T> tmpdata;

		read_netcdf_axis(filename, out_x, xname);
		read_netcdf_axis(filename, out_y, yname);
		read_netcdf_grid(filename, tmpdata, xname, yname, dataname);

		out_data.resize(tmpdata, out_y.size(), out_x.size());
		tmpdata.clear();
		return;
	}

	/**
	 * @brief      保存 Netcdf 格式的规则网格数据文件
	 *
	 * @param[in]  filename  输出文件名（无后缀）
	 * @param[in]  out_data  输出数组（注意需将网格数据保存为一个一维数据，左下到右上，行优先排序）
	 * @param[in]  xnum      网格的 x 方向顶点数量
	 * @param[in]  ynum      网格的 y 方向顶点数量
	 * @param[in]  xmin      网格的 x 方向最小值
	 * @param[in]  dx        网格的 x 方向的顶点间隔
	 * @param[in]  ymin      网格的 y 方向最小值
	 * @param[in]  dy        网格的 y 方向的顶点间隔
	 * @param[in]  zmin      网格数据的最小值，默认为 NAN 。 此时函数会自动计算数组的最小值
	 * @param[in]  zmax      网格数据的最大值，默认为 NAN 。 此时函数会自动计算数组的最大值
	 * @param[in]  xname     文件中 x 坐标轴的变量名，默认为 x
	 * @param[in]  yname     文件中 y 坐标轴的变量名，默认为 y
	 * @param[in]  dataname  文件中网格数据的变量名，默认为 z
	 */
	template <typename T, typename P>
	void save_netcdf_grid(std::string filename, const array<T> &out_data, int xnum, int ynum,
		P xmin, P dx, P ymin, P dy, std::string xname = "x", std::string yname = "y", 
		std::string dataname = "z", T zmin = NAN, T zmax = NAN)
	{
		// 准备写出数据
		if (out_data.empty())
		{
			throw runtime_error("The output data is empty. From save_netcdf_grid(...)");
		}

		if (filename == "")
		{
			throw domain_error("The input filename is empty. From save_netcdf_grid(...)");
		}

		// 判断文件名结尾是否为.nc如果不是则添加.nc结尾
		std::string full_name;
		if (filename.length() <= 3 || filename.substr(filename.length()-3, filename.length()) != ".nc")
		{
			full_name = filename + ".nc";
		}
		else full_name = filename;

		NcFile nc_outfile(full_name.c_str(), NcFile::Replace);
		if (!nc_outfile.is_valid())
		{
			throw runtime_error("Can't not open file: " + full_name + ". From save_netcdf_grid(...)");
		}

		if (std::isnan(zmin))
		{
			zmin = GCTL_BDL_MAX;
			for (int j = 0; j < out_data.size(); j++)
			{
				if (!std::isnan(out_data[j]) && out_data[j] != GCTL_BDL_MAX)
					zmin = GCTL_MIN(zmin, out_data[j]);
			}
		}

		if (std::isnan(zmax))
		{
			zmax = GCTL_BDL_MIN;
			for (int j = 0; j < out_data.size(); j++)
			{
				if (!std::isnan(out_data[j]) && out_data[j] != GCTL_BDL_MAX)
					zmax = GCTL_MAX(zmax, out_data[j]);
			}
		}

		T data_range[2] = {zmin, zmax};
		P x_range[2] = {xmin, xmin+(xnum-1)*dx};
		P y_range[2] = {ymin, ymin+(ynum-1)*dy};
		array<P> xs(xnum), ys(ynum);
		for (int i = 0; i < xnum; i++)
		{
			xs[i] = xmin + i*dx;
		}

		for (int i = 0; i < ynum; i++)
		{
			ys[i] = ymin + i*dy;
		}

		time_t now = time(0);
		char* dt_char = ctime(&now);
		dt_char[strlen(dt_char)-1] = 0; // 删除换行符
		std::string dt = dt_char;
		std::string his_str = "Generated by gctl::save_netcdf_grid() on " + dt;
		//nc_outfile.add_att("title", "Topography data.");
		nc_outfile.add_att("history", his_str.c_str());

		NcDim *xDim = nc_outfile.add_dim(xname.c_str(), xnum);
		NcDim *yDim = nc_outfile.add_dim(yname.c_str(), ynum);

		NcType axis_type, val_type;
		if (typeid(P).name() == typeid(signed char).name()) axis_type = ncByte;
		else if (typeid(P).name() == typeid(char).name()) axis_type = ncChar;
		else if (typeid(P).name() == typeid(short).name()) axis_type = ncShort;
		else if (typeid(P).name() == typeid(int).name()) axis_type = ncInt;
		else if (typeid(P).name() == typeid(float).name()) axis_type = ncFloat;
		else if (typeid(P).name() == typeid(double).name()) axis_type = ncDouble;
		else axis_type = ncDouble;

		if (typeid(T).name() == typeid(signed char).name()) val_type = ncByte;
		else if (typeid(T).name() == typeid(char).name()) val_type = ncChar;
		else if (typeid(T).name() == typeid(short).name()) val_type = ncShort;
		else if (typeid(T).name() == typeid(int).name()) val_type = ncInt;
		else if (typeid(T).name() == typeid(float).name()) val_type = ncFloat;
		else if (typeid(T).name() == typeid(double).name()) val_type = ncDouble;
		else val_type = ncDouble;

		NcVar *xVar = nc_outfile.add_var(xname.c_str(), axis_type, xDim);
		//xVar->add_att("units", 'm');
		xVar->add_att("valid_range", 2, &x_range[0]);
		xVar->put(xs.get(), xnum);

		NcVar *yVar = nc_outfile.add_var(yname.c_str(), axis_type, yDim);
		//yVar->add_att("units", 'm');
		yVar->add_att("valid_range", 2, &y_range[0]);
		yVar->put(ys.get(), ynum);

		NcVar *dVar = nc_outfile.add_var(dataname.c_str(), val_type, yDim, xDim);
		//dVar->add_att("units", 'm');
		dVar->add_att("valid_range", 2, &data_range[0]);
		dVar->add_att("missing_value", GCTL_BDL_MAX);

		array<T> tmp_data(out_data);
		for (int i = 0; i < tmp_data.size(); i++)
		{
			if (tmp_data[i] < zmin || tmp_data[i] > zmax)
				tmp_data[i] = GCTL_BDL_MAX;
		}
		dVar->put(tmp_data.get(), ynum, xnum);
		return;
	}

	/**
	 * @brief      保存 Netcdf 格式的规则网格数据文件
	 *
	 * @param[in]  filename  输出文件名（无后缀）
	 * @param[in]  out_data  输出二维数组
	 * @param[in]  xmin      网格的 x 方向最小值
	 * @param[in]  dx        网格的 x 方向的顶点间隔
	 * @param[in]  ymin      网格的 y 方向最小值
	 * @param[in]  dy        网格的 y 方向的顶点间隔
	 * @param[in]  zmin      网格数据的最小值，默认为 NAN 。 此时函数会自动计算数组的最小值
	 * @param[in]  zmax      网格数据的最大值，默认为 NAN 。 此时函数会自动计算数组的最大值
	 * @param[in]  xname     文件中 x 坐标轴的变量名，默认为 x
	 * @param[in]  yname     文件中 y 坐标轴的变量名，默认为 y
	 * @param[in]  dataname  文件中网格数据的变量名，默认为 z
	 */
	template <typename T, typename P>
	void save_netcdf_grid(std::string filename, const matrix<T> &out_data, 
		P xmin, P dx, P ymin, P dy, std::string xname = "x", std::string yname = "y", 
		std::string dataname = "z", T zmin = NAN, T zmax = NAN)
	{
		// 准备写出数据
		if (out_data.empty())
		{
			throw runtime_error("The output data is empty. From save_netcdf_grid(...)");
		}

		gctl::array<T> tmpdata;
		out_data.reform(tmpdata);

		save_netcdf_grid(filename, tmpdata, out_data.col_size(), out_data.row_size(), 
			xmin, dx, ymin, dy, xname, yname, dataname, zmin, zmax);
		return;
	}

	/**
	 * @brief      保存 Netcdf 格式的规则网格数据文件
	 *
	 * @param[in]  filename  输出文件名（无后缀）
	 * @param[in]  out_data  输出数组（注意需将网格数据保存为一个一维数据，左下到右上，行优先排序）
	 * @param[in]  out_x     输出的 x 坐标数组
	 * @param[in]  out_y     输出的 y 坐标数组
	 * @param[in]  zmin      网格数据的最小值，默认为 NAN 。 此时函数会自动计算数组的最小值
	 * @param[in]  zmax      网格数据的最大值，默认为 NAN 。 此时函数会自动计算数组的最大值
	 * @param[in]  xname     文件中 x 坐标轴的变量名，默认为 x
	 * @param[in]  yname     文件中 y 坐标轴的变量名，默认为 y
	 * @param[in]  dataname  文件中网格数据的变量名，默认为 z
	 */
	template <typename T, typename P>
	void save_netcdf_grid(std::string filename, const array<T> &out_data, const array<P> &out_x, 
		const array<P> &out_y, std::string xname = "x", std::string yname = "y", std::string dataname = "z", 
		T zmin = NAN, T zmax = NAN)
	{
		// 准备写出数据
		if (out_data.empty())
		{
			throw runtime_error("The output data is empty. From save_netcdf_grid(...)");
		}

		if (filename == "")
		{
			throw domain_error("The input filename is empty. From save_netcdf_grid(...)");
		}

		// 判断文件名结尾是否为.nc如果不是则添加.nc结尾
		std::string full_name;
		if (filename.length() <= 3 || filename.substr(filename.length()-3, filename.length()) != ".nc")
		{
			full_name = filename + ".nc";
		}
		else full_name = filename;

		NcFile nc_outfile(full_name.c_str(), NcFile::Replace);
		if (!nc_outfile.is_valid())
		{
			throw runtime_error("Can't not open file: " + full_name + ". From save_netcdf_grid(...)");
		}

		if (std::isnan(zmin))
		{
			zmin = GCTL_BDL_MAX;
			for (int j = 0; j < out_data.size(); j++)
			{
				if (!std::isnan(out_data[j]) && out_data[j] != GCTL_BDL_MAX)
					zmin = GCTL_MIN(zmin, out_data[j]);
			}
		}

		if (std::isnan(zmax))
		{
			zmax = GCTL_BDL_MIN;
			for (int j = 0; j < out_data.size(); j++)
			{
				if (!std::isnan(out_data[j]) && out_data[j] != GCTL_BDL_MAX)
					zmax = GCTL_MAX(zmax, out_data[j]);
			}
		}

		int xnum = out_x.size(), ynum = out_y.size();
		T data_range[2] = {zmin, zmax};
		P x_range[2] = {out_x[0], out_x[xnum-1]};
		P y_range[2] = {out_y[0], out_y[ynum-1]};

		time_t now = time(0);
		char* dt_char = ctime(&now);
		dt_char[strlen(dt_char)-1] = 0; // 删除换行符
		std::string dt = dt_char;
		std::string his_str = "Generated by gctl::save_netcdf_grid() on " + dt;
		//nc_outfile.add_att("title", "Topography data.");
		nc_outfile.add_att("history", his_str.c_str());

		NcDim *xDim = nc_outfile.add_dim(xname.c_str(), xnum);
		NcDim *yDim = nc_outfile.add_dim(yname.c_str(), ynum);

		NcType axis_type = ncDouble, val_type = ncDouble;
		if (typeid(P).name() == typeid(signed char).name()) axis_type = ncByte;
		else if (typeid(P).name() == typeid(char).name()) axis_type = ncChar;
		else if (typeid(P).name() == typeid(short).name()) axis_type = ncShort;
		else if (typeid(P).name() == typeid(int).name()) axis_type = ncInt;
		else if (typeid(P).name() == typeid(float).name()) axis_type = ncFloat;
		else if (typeid(P).name() == typeid(double).name()) axis_type = ncDouble;

		if (typeid(T).name() == typeid(signed char).name()) val_type = ncByte;
		else if (typeid(T).name() == typeid(char).name()) val_type = ncChar;
		else if (typeid(T).name() == typeid(short).name()) val_type = ncShort;
		else if (typeid(T).name() == typeid(int).name()) val_type = ncInt;
		else if (typeid(T).name() == typeid(float).name()) val_type = ncFloat;
		else if (typeid(T).name() == typeid(double).name()) val_type = ncDouble;

		NcVar *xVar = nc_outfile.add_var(xname.c_str(), axis_type, xDim);
		//xVar->add_att("units", 'm');
		xVar->add_att("valid_range", 2, &x_range[0]);
		xVar->put(out_x.get(), xnum);

		NcVar *yVar = nc_outfile.add_var(yname.c_str(), axis_type, yDim);
		//yVar->add_att("units", 'm');
		yVar->add_att("valid_range", 2, &y_range[0]);
		yVar->put(out_y.get(), ynum);

		NcVar *dVar = nc_outfile.add_var(dataname.c_str(), val_type, yDim, xDim);
		//dVar->add_att("units", 'm');
		dVar->add_att("valid_range", 2, &data_range[0]);
		dVar->add_att("missing_value", GCTL_BDL_MAX);

		array<T> tmp_data(out_data.get(), out_data.size());
		for (int i = 0; i < tmp_data.size(); i++)
		{
			if (tmp_data[i] < zmin || tmp_data[i] > zmax)
				tmp_data[i] = GCTL_BDL_MAX;
		}
		dVar->put(tmp_data.get(), ynum, xnum);
		return;
	}

	/**
	 * @brief      保存 Netcdf 格式的规则网格数据文件
	 *
	 * @param[in]  filename  输出文件名（无后缀）
	 * @param[in]  out_data  输出二维数组
	 * @param[in]  out_x     输出的 x 坐标数组
	 * @param[in]  out_y     输出的 y 坐标数组
	 * @param[in]  zmin      网格数据的最小值，默认为 NAN 。 此时函数会自动计算数组的最小值
	 * @param[in]  zmax      网格数据的最大值，默认为 NAN 。 此时函数会自动计算数组的最大值
	 * @param[in]  xname     文件中 x 坐标轴的变量名，默认为 x
	 * @param[in]  yname     文件中 y 坐标轴的变量名，默认为 y
	 * @param[in]  dataname  文件中网格数据的变量名，默认为 z
	 */
	template <typename T, typename P>
	void save_netcdf_grid(std::string filename, const matrix<T> &out_data, const array<P> &out_x, 
		const array<P> &out_y, std::string xname = "x", std::string yname = "y", std::string dataname = "z", 
		T zmin = NAN, T zmax = NAN)
	{
		// 准备写出数据
		if (out_data.empty())
		{
			throw runtime_error("The output data is empty. From save_netcdf_grid(...)");
		}

		gctl::array<T> tmpdata;
		out_data.reform(tmpdata);

		save_netcdf_grid(filename, tmpdata, out_x, out_y, xname, yname, dataname, zmin, zmax);
		return;
	}

	/**
	 * @brief      追加保存 Netcdf 格式的规则网格数据文件
	 * 
	 * @warning    追加的网格数据大小必须与已有网格一致
	 *
	 * @param[in]  filename  输出文件名（无后缀）
	 * @param[in]  out_data  输出数组（注意需将网格数据保存为一个一维数据，左下到右上，行优先排序）
	 * @param[in]  xname     保存到的 x 维度名称
	 * @param[in]  yname     保存到的 y 维度名称
	 * @param[in]  dataname  保存的数据名称，注意不要与已有数据重名
	 * @param[in]  zmin      网格数据的最小值，默认为 NAN 。 此时函数会自动计算数组的最小值
	 * @param[in]  zmax      网格数据的最大值，默认为 NAN 。 此时函数会自动计算数组的最大值
	 */
	template <typename T>
	void append_netcdf_grid(std::string filename, const array<T> &out_data, std::string xname, 
		std::string yname, std::string dataname, T zmin = NAN, T zmax = NAN)
	{
		// 准备写出数据
		if (out_data.empty())
		{
			throw runtime_error("The output data is empty. From append_netcdf_grid(...)");
		}

		if (filename == "")
		{
			throw domain_error("The input filename is empty. From append_netcdf_grid(...)");
		}

		// 判断文件名结尾是否为.nc如果不是则添加.nc结尾
		std::string full_name;
		if (filename.length() <= 3 || filename.substr(filename.length()-3, filename.length()) != ".nc")
		{
			full_name = filename + ".nc";
		}
		else full_name = filename;

		NcFile nc_outfile(full_name.c_str(), NcFile::Write);
		if (!nc_outfile.is_valid())
		{
			throw runtime_error("Can't not open file: " + full_name + ". From append_netcdf_grid(...)");
		}

		NcVar *tmp_var;
		std::string tmp_name;
		for (int i = 0; i < nc_outfile.num_vars(); i++)
		{
			tmp_var = nc_outfile.get_var(i);
			tmp_name = tmp_var->name();
			if (tmp_name == dataname)
			{
				throw runtime_error("Variable already existed: " + dataname + ". From append_netcdf_grid(...)");
			}
		}

		if (std::isnan(zmin))
		{
			zmin = GCTL_BDL_MAX;
			for (int j = 0; j < out_data.size(); j++)
			{
				if (!std::isnan(out_data[j]) && out_data[j] != GCTL_BDL_MAX)
					zmin = GCTL_MIN(zmin, out_data[j]);
			}
		}

		if (std::isnan(zmax))
		{
			zmax = GCTL_BDL_MIN;
			for (int j = 0; j < out_data.size(); j++)
			{
				if (!std::isnan(out_data[j]) && out_data[j] != GCTL_BDL_MAX)
					zmax = GCTL_MAX(zmax, out_data[j]);
			}
		}

		T data_range[2] = {zmin, zmax};

		time_t now = time(0);
		char* dt_char = ctime(&now);
		dt_char[strlen(dt_char)-1] = 0; // 删除换行符
		std::string dt = dt_char;
		std::string his_str = "Edited by gctl::append_netcdf_grid() on " + dt;
		//nc_outfile.add_att("title", "Topography data.");
		nc_outfile.add_att("history", his_str.c_str());

		NcDim *xDim = nc_outfile.get_dim(xname.c_str());
		NcDim *yDim = nc_outfile.get_dim(yname.c_str());

		if (!xDim->is_valid())
		{
			throw runtime_error("Dimension not found: " + xname + ". From append_netcdf_grid(...)");
		}

		if (!yDim->is_valid())
		{
			throw runtime_error("Dimension not found: " + yname + ". From append_netcdf_grid(...)");
		}

		int xnum = xDim->size();
		int ynum = yDim->size();

		if (xnum * ynum != out_data.size())
		{
			throw runtime_error("Grid size does not match. From append_netcdf_grid(...)");
		}

		NcType val_type = ncDouble;
		if (typeid(T).name() == typeid(signed char).name()) val_type = ncByte;
		else if (typeid(T).name() == typeid(char).name()) val_type = ncChar;
		else if (typeid(T).name() == typeid(short).name()) val_type = ncShort;
		else if (typeid(T).name() == typeid(int).name()) val_type = ncInt;
		else if (typeid(T).name() == typeid(float).name()) val_type = ncFloat;
		else if (typeid(T).name() == typeid(double).name()) val_type = ncDouble;

		NcVar *dVar = nc_outfile.add_var(dataname.c_str(), val_type, yDim, xDim);
		//dVar->add_att("units", 'm');
		dVar->add_att("valid_range", 2, &data_range[0]);
		dVar->add_att("missing_value", GCTL_BDL_MAX);

		array<T> tmp_data(out_data.get(), out_data.size());
		for (int i = 0; i < tmp_data.size(); i++)
		{
			if (tmp_data[i] < zmin || tmp_data[i] > zmax)
				tmp_data[i] = GCTL_BDL_MAX;
		}
		dVar->put(tmp_data.get(), ynum, xnum);

		return;
	}

	/**
	 * @brief      追加保存 Netcdf 格式的规则网格数据文件
	 * 
	 * @warning    追加的网格数据大小必须与已有网格一致
	 *
	 * @param[in]  filename  输出文件名（无后缀）
	 * @param[in]  out_data  输出二维数组
	 * @param[in]  xname     保存到的 x 维度名称
	 * @param[in]  yname     保存到的 y 维度名称
	 * @param[in]  dataname  保存的数据名称，注意不要与已有数据重名
	 * @param[in]  zmin      网格数据的最小值，默认为 NAN 。 此时函数会自动计算数组的最小值
	 * @param[in]  zmax      网格数据的最大值，默认为 NAN 。 此时函数会自动计算数组的最大值
	 */
	template <typename T>
	void append_netcdf_grid(std::string filename, const matrix<T> &out_data, std::string xname, 
		std::string yname, std::string dataname, T zmin = NAN, T zmax = NAN)
	{
		// 准备写出数据
		if (out_data.empty())
		{
			throw runtime_error("The output data is empty. From append_netcdf_grid(...)");
		}

		gctl::array<T> tmpdata;
		out_data.reform(tmpdata);

		append_netcdf_grid(filename, tmpdata, xname, yname, dataname, zmin, zmax);

		return;
	}
}

#endif // GCTL_NETCDF

#endif // _GCTL_NETCDF_IO_H