最近运用到QGis的开发,通过gdal库(操作各种栅格地理数据格式的库)解析shp文件得到想要的地理数据,开发过程中,想得到某段道路的折线长度,网上找很多资料都没有类似的代码示例(主要是讲解地理坐标和投影坐标的定义,很清晰、易理解,但没有找到自己想要的示例代码)。
开发环境:visual studio2013 (windows) gdal-2.3.1
例:需要得到图中这段道路的长度(gdal开发)
如图所示,gdal代码需要获取道路的长度921.130米
OGRLineString *lineString;
lineString->get_Length();//获取的长度为笛卡尔坐标长度0.011deg
查询发现,发现好像没有可以直接获取椭球体的长度接口,因此获取的点坐标需要进行转换,转换代码很简单。
QGis软件所显示epsg:4326是地理坐标,需要的坐标为投影坐标
(通过postgresql数据库开发的话,可以直接添加命令ST_Length(ST_TransFrom(geom,4326),true)获得)
转换代码:
//points参数为获取的折线的所有点坐标,存放在一个vector中
double GdalDevelop::calculateLength(vector<OGRPoint> points)
{
char *text = nullptr;
OGRSpatialReference source, target;
OGRErr err1 = source.importFromEPSG(4326);
OGRErr err2 = target.importFromEPSG(32650);//32650通过qgis软件查看得到epsg
source.exportToWkt(&text);//导出打印wkt
cout <<"text1:" <<text << endl;
target.exportToWkt(&text);//导出打印wkt
cout <<"text2:" <<text<<endl;
if (err1 != OGRERR_NONE || err2 != OGRERR_NONE)
{
cout << err1 << " "<<err2 << endl;
return 0;
}
//初始化转化类
OGRCoordinateTransformation *coord = OGRCreateCoordinateTransformation(&source, &target);
//根据输入的点计算折线长度
double total_len = 0;
int num = points.size();
for (int i = 0; i < num; i++)
{
if (i+1 < num)
{
double x0 = points[i].getX();
double y0 = points[i].getY();
double x1 = points[i + 1].getX();
double y1 = points[i + 1].getY();
//cout << "转换前:" << x0 << " " << y0 << endl;
coord->Transform(1, &x0, &y0);
//cout << "转换后:" << x0 << " " << y0 << endl;
//cout << "转换前:" << x1 << " " << y1 << endl;
coord->Transform(1, &x1, &y1);
//cout << "转换后:" << x1 << " " << y1 << endl;
double diff_x = x1 - x0;
double diff_y = y1 - y0;
//cout <<"diff:" <<diff_x <<" "<< diff_y<<endl;
total_len += sqrt(diff_x*diff_x + diff_y*diff_y);
}
}
return total_len;
}
//计算得到的结果是存在点点误差的
32650的由来:
需要得到投影坐标对应的坐标参考系如图:WGS 84/UTM zone 50N, 选择坐标系过程中需要一一对应,如果不对图层覆盖范围是不对的。
GDAL投影转换:https://www.cnblogs.com/renyu310/p/5907072.html
QGIS视频网址:https://www.bilibili.com/video/BV1op4y1D7PS?from=search&seid=14366299219202295850