I can not deal with a strange bug. At the input: the coordinate of the WGS84 point and the radius in meters. The task is quite simple, you need to build a buffer around this point so that it is correctly displayed in WGS84 coordinates.
Solution:
- create an Azimuthal Equidistant projection with the center at the specified point
- using the GDAL library and the standard "buffer" function, create a polygon
- reproject the polygon back to WGS84
Where is the mistake:
- coordinates are truncated
- there is an incorrect transformation of the polygon
For reprojection, the PROJ library and GeographicLib::AzimuthalEquidistant were used.
(+proj=aeqd +lon_0=129.4943 +lat_0=64.2402 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs +type=crs
)
The result is the same.
Moreover, if the task is completed in ArcGIS, then the ideal option is obtained.
Code example:
void convert_wgs84_to_aeqd(double lonCentr, double latCentr, double lonInp, double latInp, double& lonRes, double& latRes, bool isToWgs)
{
Geodesic geod(Constants::WGS84_a(), Constants::WGS84_f());
AzimuthalEquidistant proj(geod);
const double lat0 = latCentr, lon0 = lonCentr;
if (!isToWgs) {
proj.Forward(lat0, lon0, latInp, lonInp, lonRes, latRes);
}
else {
proj.Reverse(lat0, lon0, latInp, lonInp, latRes, lonRes);
}
}
void create_buffer(double lonCentr, double latCentr, double bufferInMeters)
{
const double buffer_distance = bufferInMeters;
const int points_per_circle = 360;
using coordinate_type = double;
boost::geometry::strategy::buffer::distance_symmetric<coordinate_type> distance_strategy(buffer_distance);
boost::geometry::strategy::buffer::join_round join_strategy(points_per_circle);
boost::geometry::strategy::buffer::end_round end_strategy(points_per_circle);
boost::geometry::strategy::buffer::point_circle circle_strategy(points_per_circle);
boost::geometry::strategy::buffer::side_straight side_strategy;
double xr, yr;
convert_wgs84_to_aeqd(lonCentr, latCentr, lonCentr, latCentr, xr, yr, false);
boost_point px(xr, yr);
boost_mulpolygon polyx;
bg::buffer(px, polyx, distance_strategy, side_strategy, join_strategy, end_strategy, circle_strategy);
OGRLinearRing* glr = new OGRLinearRing;
auto rng = bg::exterior_ring(polyx[0]);
for (auto vp : rng) {
double laty = bg::get<1>(vp);
double lonx = bg::get<0>(vp);
double lonRes, latRes;
convert_wgs84_to_aeqd(lonCentr, latCentr, lonx, laty, lonRes, latRes, true);
glr->addPoint(lonRes, latRes);
}
OGRPolygon* pdff = new OGRPolygon;
pdff->addRingDirectly(glr);
char* wktStrings;
pdff->exportToWkt(&wktStrings);
std::string r(wktStrings);
std::cout<< r <<std::endl;
}