(GDAL-OGR)利用EPSG_Code创建wkt格式的坐标系描述
上一篇中讲到如何提取las文件中的地理坐标系信息,提及EPSG_Code可以通过GDAL进行坐标系描述和转换,这里进行代码测试和详细介绍:
GDAL:
再GDAL中,创建带有地理坐标的tif栅格文件或者shp矢量文件,都需要制定坐标系,我习惯使用wkt格式坐标系描述方法,(就是将坐标系信息用wkt格式表示)
GDALDataset* pCreateDataset = pDriver->Create(strFileFullPath.c_str(),nlmgWidth,nlmgHeight,ChannelNum,dataTypeForCreate,papszOptions);
pCreateDataset->SetProjection (strWkt)
一般情况下我们可以通过如下方法获取指定坐标系的wkt描述:
char *strWkt = NULL;
if (Wgs84ORUTM == 0 )
{
OGRSpatialReference oSRS;
oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&strWkt);
int a = oSRS.IsGeographic();
int b = oSRS.IsProjected();
}
else if(Wgs84ORUTM == 1)
{
OGRSpatialReference oSRS;
oSRS.SetProjCS("UTM (WGS84) in northernhemisphere.");
oSRS.SetWellKnownGeogCS("WGS84");
oSRS.SetUTM(UTMNum,TRUE);
oSRS.exportToWkt(&strWkt);
//strWkt ="PROJCS["UTM",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG",6326]],PRIMEM["Greenwich",0,AUTHORITY["EPSG",8901]],UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG",9108]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG",4326]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0]]";
}
但是上面的定义方法比较局限,默认GCS就是WGS84,PCS就是UTM;当我们从las文件中提取地理坐标系信息时,有可能会遇到各种各样的坐标系,为了提高接口的兼容性,我们利用EPSG_Code来传递坐标系信息,我们将las文件中获取的EPSG_Code传递给GDAL,然后用下面的方法获取wkt格式的坐标系描述,就可以生成相同坐地理坐标系的tif/shp等文件了;
OGRSpatialReference oSRS;
oSRS.importFromEPSG((int)4214);//北京54
char* strWkt;
oSRS.exportToWkt(&strWkt);
printf("%s", strWkt);
对于EPSG_Code和GDAL的说明:
通过测试和阅读说明发现,地理坐标系(GCS)最好只用GCS_开头的(既包含椭球面又指定了基准面),GCSE_开头的头没有指定基准面(Datum)只包含椭球面,测试只要GCS_开头的就能正确导入到OGR中;投影坐标系(PCS)就比较简单了,就是对应的那些code;
另外,在使用GDAL导入EPSG_Code时需要将GDAL库目录下…/data目录中的文件全部拷贝到工程目录下,否则查询不到相应的EPSG_Code:
转载自:https://blog.csdn.net/xujie126/article/details/80988850