Creating a buffer in meters, Azimuthal Equidistant projection

174 Views Asked by At

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:

  1. create an Azimuthal Equidistant projection with the center at the specified point
  2. using the GDAL library and the standard "buffer" function, create a polygon
  3. reproject the polygon back to WGS84

Where is the mistake:

  1. coordinates are truncated
  2. 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.

enter image description here enter image description here

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;
}
0

There are 0 best solutions below