235 lines
7.6 KiB
C++
235 lines
7.6 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2022 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 "gctl/core.h"
|
|
#include "gctl/io.h"
|
|
#include "gctl/utility.h"
|
|
|
|
int main(int argc, char *argv[]) try
|
|
{
|
|
gctl::flags_parser fp;
|
|
fp.set_proname("sph2plane");
|
|
fp.set_proinfo("Project spherical data points to a plane surface with respected to \
|
|
a reference point. This program is a toolkit of the GCTL package. The GCTL comes with ABSOLUTE NO WARRANTY. \
|
|
Please see instructions or contact the author for more information. The program \
|
|
takes inputs from the standard input and display the results on screen. The coordinates \
|
|
should be separated by spaces and different points should be placed on different lines. \
|
|
The inputting continues until encounter a blank line or 'EOF'. Any lines start with '#' will be skipped.");
|
|
|
|
fp.add_opt('a', "anchor", required_argument, NULL,
|
|
"The reference point (anchor) for the project.", "<lon>,<lat>", true);
|
|
fp.add_opt('r', "reference", required_argument, NULL,
|
|
"Set a reference system for converting the coordinates according to the specified spherical or ellipsoidal surface.", "<r>/<R>|<f>,<R>|WGS84|Earth|Moon", true);
|
|
fp.add_opt('u', "unit-length", no_argument, NULL,
|
|
"Use the unit-length at the reference point to project all data points. This is often required to get unified intervals for grid data points.", 0, false);
|
|
fp.add_opt('c', "column", required_argument, NULL,
|
|
"Select data columns for the longitudinal and latitudinal coordinates (the default is 0,1).", "<lon-col>,<lat-col>", false);
|
|
fp.add_opt('p', "precision", required_argument, NULL,
|
|
"Set the precision of the outputs (the default is 12).", "<precision>", false);
|
|
fp.add_opt('h', "help", no_argument, NULL, "Show help information.", 0, false);
|
|
fp.configure(argc, argv);
|
|
|
|
if(argc == 1 || fp.set_opt('h'))
|
|
{
|
|
fp.show_help_page();
|
|
return 0;
|
|
}
|
|
|
|
// 查看是否通过强制参数检查
|
|
if (!fp.pass_mandatory()) return 0;
|
|
|
|
bool use_unit = false;
|
|
if (fp.set_opt('u'))
|
|
{
|
|
use_unit = true;
|
|
}
|
|
|
|
std::string err_str, anchor_str, ref_str, col_str, pre_str;
|
|
std::stringstream tmp_ss;
|
|
|
|
fp.get_argv({'a', 'r', 'c', 'p'}, {&anchor_str, &ref_str, &col_str, &pre_str});
|
|
// 确定参考点
|
|
double ref_lon = -200.0, ref_lat = -100.0;
|
|
if (2 != sscanf(anchor_str.c_str(), "%lf,%lf", &ref_lon, &ref_lat) ||
|
|
ref_lon < -180.0 || ref_lon > 180.0 || ref_lat < -90.0 || ref_lat > 90.0)
|
|
{
|
|
err_str = "Wrong entrance: " + anchor_str;
|
|
throw err_str;
|
|
}
|
|
|
|
// 确定参考球
|
|
double ref_r, ref_R;
|
|
if (ref_str == "WGS84")
|
|
{
|
|
ref_r = GCTL_WGS84_PoleRadius;
|
|
ref_R = GCTL_WGS84_EquatorRadius;
|
|
}
|
|
else if (ref_str == "Earth")
|
|
{
|
|
ref_r = GCTL_Earth_Radius;
|
|
ref_R = GCTL_Earth_Radius;
|
|
}
|
|
else if (ref_str == "Moon")
|
|
{
|
|
ref_r = GCTL_Moon_Radius;
|
|
ref_R = GCTL_Moon_Radius;
|
|
}
|
|
else if (2 != sscanf(ref_str.c_str(), "%lf/%lf", &ref_r, &ref_R))
|
|
{
|
|
if (2 == sscanf(ref_str.c_str(), "%lf,%lf", &ref_r, &ref_R))
|
|
{
|
|
ref_r = (1.0-ref_r)*ref_R;
|
|
}
|
|
else
|
|
{
|
|
err_str = "Wrong entrance: " + ref_str;
|
|
throw err_str;
|
|
}
|
|
}
|
|
|
|
// 确定坐标列
|
|
int col_lon, col_lat;
|
|
if (col_str == "NULL") col_str = "0,1";
|
|
if (2 != sscanf(col_str.c_str(), "%d,%d", &col_lon, &col_lat))
|
|
{
|
|
err_str = "Wrong entrance: " + col_str;
|
|
throw err_str;
|
|
}
|
|
|
|
// 确定输出精度
|
|
int pre_int;
|
|
if (pre_str == "NULL") pre_str = "12";
|
|
if (1 != sscanf(pre_str.c_str(), "%d", &pre_int))
|
|
{
|
|
err_str = "Wrong entrance: " + pre_str;
|
|
throw err_str;
|
|
}
|
|
|
|
std::string line;
|
|
std::vector<std::string> line_arr;
|
|
std::vector<double> line_ele;
|
|
gctl::point3dc pc;
|
|
gctl::point3ds ps;
|
|
while (true)
|
|
{
|
|
std::getline(std::cin, line);
|
|
if (line[0] == '#') continue;
|
|
if (line.empty() || line == "EOF") break;
|
|
else line_arr.push_back(line);
|
|
}
|
|
|
|
std::cout << "# Outputs from the sph2plane. Use -h to see the help information." << std::endl;
|
|
std::cout << "# reference system: " << ref_str << std::endl;
|
|
std::cout << "# reference coordinate: lon = " << ref_lon << ", lat = " << ref_lat << std::endl;
|
|
|
|
double rad, rad_lat, unit_yl, unit_xl, tmp_d;
|
|
double point_lon, point_lat;
|
|
if (use_unit)
|
|
{
|
|
// 计算参考点处的经纬向单位长度
|
|
rad = gctl::ellipse_radius_2d(ref_R, ref_r, ref_lat*GCTL_Pi/180.0);
|
|
unit_yl = rad*tan(1.0*GCTL_Pi/180.0);
|
|
|
|
rad_lat = rad*sin(GCTL_Pi*(90.0 - ref_lat)/180.0);
|
|
unit_xl = rad_lat*tan(1.0*GCTL_Pi/180.0);
|
|
|
|
for (int i = 0; i < line_arr.size(); i++)
|
|
{
|
|
if (!line_ele.empty()) line_ele.clear();
|
|
|
|
tmp_ss.clear();
|
|
tmp_ss.str(line_arr[i]);
|
|
while (tmp_ss.good())
|
|
{
|
|
tmp_ss >> tmp_d;
|
|
line_ele.push_back(tmp_d);
|
|
}
|
|
|
|
if (line_ele.size() <= GCTL_MAX(col_lon, col_lat))
|
|
{
|
|
err_str = "Wrong entrance: " + line_arr[i];
|
|
throw err_str;
|
|
}
|
|
|
|
point_lon = line_ele[col_lon];
|
|
point_lat = line_ele[col_lat];
|
|
|
|
line_ele[col_lon] = (point_lon - ref_lon)*unit_xl;
|
|
line_ele[col_lat] = (point_lat - ref_lat)*unit_yl;
|
|
|
|
std::cout << std::setprecision(pre_int) << line_ele[0];
|
|
for (int l = 1; l < line_ele.size(); l++)
|
|
{
|
|
std::cout << std::setprecision(pre_int) << " " << line_ele[l];
|
|
}
|
|
std::cout << std::endl;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for (int i = 0; i < line_arr.size(); i++)
|
|
{
|
|
if (!line_ele.empty()) line_ele.clear();
|
|
|
|
tmp_ss.clear();
|
|
tmp_ss.str(line_arr[i]);
|
|
while (tmp_ss.good())
|
|
{
|
|
tmp_ss >> tmp_d;
|
|
line_ele.push_back(tmp_d);
|
|
}
|
|
|
|
if (line_ele.size() <= GCTL_MAX(col_lon, col_lat))
|
|
{
|
|
err_str = "Wrong entrance: " + line_arr[i];
|
|
throw err_str;
|
|
}
|
|
|
|
point_lon = line_ele[col_lon];
|
|
point_lat = line_ele[col_lat];
|
|
|
|
// 计算参考点处的经纬向单位长度
|
|
rad = gctl::ellipse_radius_2d(ref_R, ref_r, ref_lat*GCTL_Pi/180.0);
|
|
rad_lat = rad*sin(GCTL_Pi*(90.0 - ref_lat)/180.0);
|
|
|
|
line_ele[col_lat] = rad*tan((point_lat - ref_lat)*GCTL_Pi/180.0);
|
|
line_ele[col_lon] = rad_lat*tan((point_lon - ref_lon)*GCTL_Pi/180.0);
|
|
|
|
std::cout << std::setprecision(pre_int) << line_ele[0];
|
|
for (int l = 1; l < line_ele.size(); l++)
|
|
{
|
|
std::cout << std::setprecision(pre_int) << " " << line_ele[l];
|
|
}
|
|
std::cout << std::endl;
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
catch(std::exception &e)
|
|
{
|
|
GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
|
|
} |