基于Azimuthal Equidistant投影创建米级缓冲区的异常问题

输入为WGS84坐标点及米级半径，需围绕该点创建缓冲区并在WGS84坐标下正确显示，但遇到异常问题。

解决方案思路 #

• 以指定点为中心创建等距方位（Azimuthal Equidistant）投影
• 使用GDAL库及标准buffer函数生成多边形
• 将多边形重投影回WGS84

出现的错误 #

• 坐标被截断
• 多边形转换结果不正确

重投影参数 #

重投影时使用了PROJ库及GeographicLib::AzimuthalEquidistant，参数为：

+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

两种方式结果一致，但在ArcGIS中可得到理想结果。

代码示例 #

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

内容的提问来源于stack exchange，提问作者luzhskij