gctl/lib/io/surfer_io.cpp
2024-09-10 15:45:07 +08:00

660 lines
18 KiB
C++

/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* 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.
******************************************************/
#include "surfer_io.h"
void gctl::read_surfer6_grid(std::string filename, _2d_matrix &out_data,
int &xnum, int &ynum, double &xmin, double &xmax, double &ymin, double &ymax,
double &zmin, double &zmax, surfer_file_type_e file_format)
{
std::ifstream infile;
if (file_format == Surfer6Binary)
{
open_infile(infile, filename, ".grd", std::ios::in|std::ios::binary);
char sym_id[5] = "NULL";
for (int i = 0; i < 4; i++)
{
infile.read((char*)&sym_id[i], sizeof(char));
}
if (strcmp(sym_id, "DSBB"))
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer6_grid(...)");
}
float tmp_d;
short int M, N;
infile.read((char*)&N, sizeof(short int));
infile.read((char*)&M, sizeof(short int));
xnum = N;
ynum = M;
// resize space for output array
out_data.resize(M, N, 0.0);
infile.read((char*)&xmin, sizeof(double));
infile.read((char*)&xmax, sizeof(double));
infile.read((char*)&ymin, sizeof(double));
infile.read((char*)&ymax, sizeof(double));
infile.read((char*)&zmin, sizeof(double));
infile.read((char*)&zmax, sizeof(double));
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
infile.read((char*)&tmp_d, sizeof(float));
out_data[i][j] = (double)tmp_d;
}
}
infile.close();
}
else
{
open_infile(infile, filename, ".grd");
std::string tmp_str;
// read first line
infile >> tmp_str;
if (tmp_str != "DSAA")
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer6_grid(...)");
}
int M, N; // N for x number, M for y number
// read dimensions
infile >> N >> M;
xnum = N;
ynum = M;
// resize space for output array
out_data.resize(M, N, 0.0);
// read grid's dimension
infile >> xmin >> xmax;
infile >> ymin >> ymax;
infile >> zmin >> zmax;
// read grid's data
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
infile >> tmp_str;
out_data[i][j] = str2double(tmp_str);
}
}
infile.close();
}
return;
}
void gctl::read_surfer6_grid(std::string filename, array<double> &out_data, int &xnum,
int &ynum, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin,
double &zmax, surfer_file_type_e file_format)
{
std::ifstream infile;
if (file_format == Surfer6Binary)
{
open_infile(infile, filename, ".grd", std::ios::in|std::ios::binary);
char sym_id[5] = "NULL";
for (int i = 0; i < 4; i++)
{
infile.read((char*)&sym_id[i], sizeof(char));
}
if (strcmp(sym_id, "DSBB"))
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer6_grid(...)");
}
float tmp_d;
short int M, N;
infile.read((char*)&N, sizeof(short int));
infile.read((char*)&M, sizeof(short int));
xnum = N;
ynum = M;
// resize space for output array
out_data.resize(M*N, 0.0);
infile.read((char*)&xmin, sizeof(double));
infile.read((char*)&xmax, sizeof(double));
infile.read((char*)&ymin, sizeof(double));
infile.read((char*)&ymax, sizeof(double));
infile.read((char*)&zmin, sizeof(double));
infile.read((char*)&zmax, sizeof(double));
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
infile.read((char*)&tmp_d, sizeof(float));
out_data[N*i + j] = (double)tmp_d;
}
}
infile.close();
}
else
{
open_infile(infile, filename, ".grd");
std::string tmp_str;
// read first line
infile >> tmp_str;
if (tmp_str != "DSAA")
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer6_grid(...)");
}
int M, N; // N for x number, M for y number
// read dimensions
infile >> N >> M;
xnum = N;
ynum = M;
// resize space for output array
out_data.resize(M*N, 0.0);
// read grid's dimension
infile >> xmin >> xmax;
infile >> ymin >> ymax;
infile >> zmin >> zmax;
// read grid's data
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
infile >> tmp_str;
out_data[N*i + j] = str2double(tmp_str);
}
}
infile.close();
}
return;
}
void gctl::save_surfer6_grid(std::string filename, const _2d_matrix &out_data,
double xmin, double xmax, double ymin, double ymax, double zmin, double zmax,
surfer_file_type_e file_format)
{
// 准备写出数据
if (out_data.empty())
{
throw runtime_error("The output data is empty. From save_surfer6_grid(...)");
}
if (std::isnan(zmin))
{
zmin = GCTL_BDL_MAX;
for (int i = 0; i < out_data.row_size(); i++)
{
for (int j = 0; j < out_data.col_size(); j++)
{
zmin = GCTL_MIN(zmin, out_data[i][j]);
}
}
}
if (std::isnan(zmax))
{
zmax = GCTL_BDL_MIN;
for (int i = 0; i < out_data.row_size(); i++)
{
for (int j = 0; j < out_data.col_size(); j++)
{
zmax = GCTL_MAX(zmax, out_data[i][j]);
}
}
}
short int M = out_data.row_size(); // 占两个字节
short int N = out_data.col_size();
std::ofstream outfile;
if (file_format == Surfer6Binary)
{
float tmp_d;
char sym_id[] = "DSBB";
open_outfile(outfile, filename, ".grd", std::ios::out|std::ios::binary);
outfile.write((char*)&sym_id[0], sizeof(char)); // 注意不要直接输出sym_id 这里只能输出4个字节
outfile.write((char*)&sym_id[1], sizeof(char));
outfile.write((char*)&sym_id[2], sizeof(char));
outfile.write((char*)&sym_id[3], sizeof(char));
outfile.write((char*)&N, sizeof(short int));
outfile.write((char*)&M, sizeof(short int));
outfile.write((char*)&xmin, sizeof(double)); // 占8个字节
outfile.write((char*)&xmax, sizeof(double));
outfile.write((char*)&ymin, sizeof(double));
outfile.write((char*)&ymax, sizeof(double));
outfile.write((char*)&zmin, sizeof(double));
outfile.write((char*)&zmax, sizeof(double));
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
tmp_d = (float)out_data[i][j]; // 注意需要转换为float类型 占两个字节
outfile.write((char*)&tmp_d, sizeof(float));
}
}
outfile.close();
}
else
{
open_outfile(outfile, filename, ".grd");
outfile << "DSAA" << std::endl;
outfile << N << " " << M << std::endl;
outfile << xmin << " " << xmax << std::endl;
outfile << ymin << " " << ymax << std::endl;
outfile << zmin << " " << zmax << std::endl;
for (int i = 0; i < M; i++)
{
if (zmin <= out_data[i][0] && out_data[i][0] <= zmax)
outfile << out_data[i][0];
else outfile << "nan";
for (int j = 1; j < N; j++)
{
if (zmin <= out_data[i][j] && out_data[i][j] <= zmax)
outfile << " " << out_data[i][j];
else outfile << " nan";
}
outfile << std::endl;
}
outfile.close();
}
return;
}
void gctl::save_surfer6_grid(std::string filename, const array<double> &out_data,
int xnum, int ynum, double xmin, double xmax, double ymin, double ymax, double zmin,
double zmax, surfer_file_type_e file_format)
{
// 准备写出数据
if (out_data.empty())
{
throw runtime_error("The output data is empty. From save_surfer6_grid(...)");
}
if (std::isnan(zmin))
{
zmin = GCTL_BDL_MAX;
for (int i = 0; i < out_data.size(); i++)
{
zmin = GCTL_MIN(zmin, out_data[i]);
}
}
if (std::isnan(zmax))
{
zmax = GCTL_BDL_MIN;
for (int i = 0; i < out_data.size(); i++)
{
zmax = GCTL_MAX(zmax, out_data[i]);
}
}
short int M = ynum; // 占两个字节
short int N = xnum;
std::ofstream outfile;
if (file_format == Surfer6Binary)
{
float tmp_d;
char sym_id[] = "DSBB";
open_outfile(outfile, filename, ".grd", std::ios::out|std::ios::binary);
outfile.write((char*)&sym_id[0], sizeof(char)); // 注意不要直接输出sym_id 这里只能输出4个字节
outfile.write((char*)&sym_id[1], sizeof(char));
outfile.write((char*)&sym_id[2], sizeof(char));
outfile.write((char*)&sym_id[3], sizeof(char));
outfile.write((char*)&N, sizeof(short int));
outfile.write((char*)&M, sizeof(short int));
outfile.write((char*)&xmin, sizeof(double)); // 占8个字节
outfile.write((char*)&xmax, sizeof(double));
outfile.write((char*)&ymin, sizeof(double));
outfile.write((char*)&ymax, sizeof(double));
outfile.write((char*)&zmin, sizeof(double));
outfile.write((char*)&zmax, sizeof(double));
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
tmp_d = (float)out_data[N*i+j]; // 注意需要转换为float类型 占两个字节
outfile.write((char*)&tmp_d, sizeof(float));
}
}
outfile.close();
}
else
{
open_outfile(outfile, filename, ".grd");
outfile << "DSAA" << std::endl;
outfile << N << " " << M << std::endl;
outfile << xmin << " " << xmax << std::endl;
outfile << ymin << " " << ymax << std::endl;
outfile << zmin << " " << zmax << std::endl;
for (int i = 0; i < M; i++)
{
if (zmin <= out_data[N*i] && out_data[N*i] <= zmax)
outfile << out_data[N*i];
else outfile << "nan";
for (int j = 1; j < N; j++)
{
if (zmin <= out_data[N*i+j] && out_data[N*i+j] <= zmax)
outfile << " " << out_data[N*i+j];
else outfile << " nan";
}
outfile << std::endl;
}
outfile.close();
}
return;
}
void gctl::read_surfer7_grid(std::string filename, _2d_matrix &out_data, double &xmin,
double &ymin, double &dx, double &dy, double &zmin, double &zmax, double &blank_val)
{
std::ifstream infile;
open_infile(infile, filename, ".grd", std::ios::in|std::ios::binary);
int head_id, head_sec_size, head_version;
infile.read((char*)&head_id, sizeof(int));
infile.read((char*)&head_sec_size, sizeof(int));
infile.read((char*)&head_version, sizeof(int));
if (head_id != 0x42525344 || head_sec_size != 4)
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
}
int grid_id, grid_sec_size, M, N;
double rotation;
infile.read((char*)&grid_id, sizeof(int));
infile.read((char*)&grid_sec_size, sizeof(int));
infile.read((char*)&M, sizeof(int));
infile.read((char*)&N, sizeof(int));
if (grid_id != 0x44495247 || grid_sec_size != 72)
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
}
// resize space for output array
out_data.resize(M, N, 0.0);
infile.read((char*)&xmin, sizeof(double));
infile.read((char*)&ymin, sizeof(double));
infile.read((char*)&dx, sizeof(double));
infile.read((char*)&dy, sizeof(double));
infile.read((char*)&zmin, sizeof(double));
infile.read((char*)&zmax, sizeof(double));
infile.read((char*)&rotation, sizeof(double)); // 暂时无用
infile.read((char*)&blank_val, sizeof(double));
int data_id, data_sec_size;
infile.read((char*)&data_id, sizeof(int));
infile.read((char*)&data_sec_size, sizeof(int));
if (data_id != 0x41544144 || data_sec_size != M*N*8)
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
}
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
infile.read((char*)&out_data[i][j], sizeof(double));
}
}
infile.close();
return;
}
void gctl::read_surfer7_grid(std::string filename, array<double> &out_data, int &xnum,
int &ynum, double &xmin, double &ymin, double &dx, double &dy, double &zmin,
double &zmax, double &blank_val)
{
std::ifstream infile;
open_infile(infile, filename, ".grd", std::ios::in|std::ios::binary);
int head_id, head_sec_size, head_version;
infile.read((char*)&head_id, sizeof(int));
infile.read((char*)&head_sec_size, sizeof(int));
infile.read((char*)&head_version, sizeof(int));
if (head_id != 0x42525344 || head_sec_size != 4)
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
}
int grid_id, grid_sec_size, M, N;
double rotation;
infile.read((char*)&grid_id, sizeof(int));
infile.read((char*)&grid_sec_size, sizeof(int));
infile.read((char*)&M, sizeof(int));
infile.read((char*)&N, sizeof(int));
if (grid_id != 0x44495247 || grid_sec_size != 72)
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
}
xnum = N;
ynum = M;
// resize space for output array
out_data.resize(M*N, 0.0);
infile.read((char*)&xmin, sizeof(double));
infile.read((char*)&ymin, sizeof(double));
infile.read((char*)&dx, sizeof(double));
infile.read((char*)&dy, sizeof(double));
infile.read((char*)&zmin, sizeof(double));
infile.read((char*)&zmax, sizeof(double));
infile.read((char*)&rotation, sizeof(double)); // 暂时无用
infile.read((char*)&blank_val, sizeof(double));
int data_id, data_sec_size;
infile.read((char*)&data_id, sizeof(int));
infile.read((char*)&data_sec_size, sizeof(int));
if (data_id != 0x41544144 || data_sec_size != M*N*8)
{
infile.close();
throw runtime_error("Wrong file format. From read_surfer7_grid(...)");
}
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
infile.read((char*)&out_data[N*i+j], sizeof(double));
}
}
infile.close();
return;
}
void gctl::save_surfer7_grid(std::string filename, const _2d_matrix &out_data, double xmin,
double ymin, double dx, double dy, double zmin, double zmax, double blank_val)
{
// 准备写出数据
if (out_data.empty())
{
throw runtime_error("The output data is empty. From save_surfer7_grid(...)");
}
if (std::isnan(zmin))
{
zmin = GCTL_BDL_MAX;
for (int i = 0; i < out_data.row_size(); i++)
{
for (int j = 0; j < out_data.col_size(); j++)
{
zmin = GCTL_MIN(zmin, out_data[i][j]);
}
}
}
if (std::isnan(zmax))
{
zmax = GCTL_BDL_MIN;
for (int i = 0; i < out_data.row_size(); i++)
{
for (int j = 0; j < out_data.col_size(); j++)
{
zmax = GCTL_MAX(zmax, out_data[i][j]);
}
}
}
std::ofstream outfile;
open_outfile(outfile, filename, ".grd", std::ios::out|std::ios::binary);
int head_id = 0x42525344; // 固定值 查看官方说明文档
int head_sec_size = sizeof(int);
int head_version = 1;
outfile.write((char*)&head_id, sizeof(int));
outfile.write((char*)&head_sec_size, sizeof(int));
outfile.write((char*)&head_version, sizeof(int));
int grid_id = 0x44495247;
int grid_sec_size = 8*sizeof(double) + 2*sizeof(int);
int M = out_data.row_size();
int N = out_data.col_size();
double rotation = 0.0;
outfile.write((char*)&grid_id, sizeof(int));
outfile.write((char*)&grid_sec_size, sizeof(int));
outfile.write((char*)&M, sizeof(int));
outfile.write((char*)&N, sizeof(int));
outfile.write((char*)&xmin, sizeof(double));
outfile.write((char*)&ymin, sizeof(double));
outfile.write((char*)&dx, sizeof(double));
outfile.write((char*)&dy, sizeof(double));
outfile.write((char*)&zmin, sizeof(double));
outfile.write((char*)&zmax, sizeof(double));
outfile.write((char*)&rotation, sizeof(double));
outfile.write((char*)&blank_val, sizeof(double));
int data_id = 0x41544144;
int data_sec_size = M * N * 8;
double tmp_d;
outfile.write((char*)&data_id, sizeof(int));
outfile.write((char*)&data_sec_size, sizeof(int));
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
tmp_d = out_data[i][j];
outfile.write((char*)&tmp_d, sizeof(double));
}
}
outfile.close();
return;
}
void gctl::save_surfer7_grid(std::string filename, const array<double> &out_data,
int xnum, int ynum, double xmin, double ymin, double dx, double dy,
double zmin, double zmax, double blank_val)
{
// 准备写出数据
if (out_data.empty())
{
throw runtime_error("The output data is empty. From save_surfer7_grid(...)");
}
if (std::isnan(zmin))
{
zmin = GCTL_BDL_MAX;
for (int i = 0; i < out_data.size(); i++)
{
zmin = GCTL_MIN(zmin, out_data[i]);
}
}
if (std::isnan(zmax))
{
zmax = GCTL_BDL_MIN;
for (int i = 0; i < out_data.size(); i++)
{
zmax = GCTL_MAX(zmax, out_data[i]);
}
}
std::ofstream outfile;
open_outfile(outfile, filename, ".grd", std::ios::out|std::ios::binary);
int head_id = 0x42525344; // 固定值 查看官方说明文档
int head_sec_size = sizeof(int);
int head_version = 1;
outfile.write((char*)&head_id, sizeof(int));
outfile.write((char*)&head_sec_size, sizeof(int));
outfile.write((char*)&head_version, sizeof(int));
int grid_id = 0x44495247;
int grid_sec_size = 8*sizeof(double) + 2*sizeof(int);
int M = ynum;
int N = xnum;
double rotation = 0.0;
outfile.write((char*)&grid_id, sizeof(int));
outfile.write((char*)&grid_sec_size, sizeof(int));
outfile.write((char*)&M, sizeof(int));
outfile.write((char*)&N, sizeof(int));
outfile.write((char*)&xmin, sizeof(double));
outfile.write((char*)&ymin, sizeof(double));
outfile.write((char*)&dx, sizeof(double));
outfile.write((char*)&dy, sizeof(double));
outfile.write((char*)&zmin, sizeof(double));
outfile.write((char*)&zmax, sizeof(double));
outfile.write((char*)&rotation, sizeof(double));
outfile.write((char*)&blank_val, sizeof(double));
int data_id = 0x41544144;
int data_sec_size = M * N * 8;
double tmp_d;
outfile.write((char*)&data_id, sizeof(int));
outfile.write((char*)&data_sec_size, sizeof(int));
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
{
tmp_d = out_data[N*i+j];
outfile.write((char*)&tmp_d, sizeof(double));
}
}
outfile.close();
return;
}