Wechat: yu389741| Email: gisdqy@163.com

Shop:https://www.giserdqy.com/shop

GDAL创建带有空间参考坐标系的shp矢量文件遇到的问题与解决方法


      之前项目需求:创建矢量文件;

    网上查找GDAL/OGR创建矢量文件的各种博客,也查了GDAL官网手册,可是实际生成的shp总是没有成功写入坐标系,空间参考坐标系总是unknown,这就让人很困惑了,实际上困惑了好几个月,断断续续的抽时间解决这个说大不大的问题,今天终于解决了,趁热跟大家分享。

    之前的代码:

/**********************************key_code***********************************************/

        //原始栅格影响空间参考信息直接赋值给矢量文件

        const char *WKT = pInDataset->GetProjectionRef();	

	if (oSRS.SetFromUserInput(WKT ) != OGRERR_NONE)
	    {	       
		    GDALDestroyDriverManager();
	            return 6;//矢量文件写入参考空间信息失败
	    }		
	poLayer = poDS->CreateLayer((const char*)FireArea_SHPpath, InputOGRS, wkbUnknown, NULL);

/*************************************key_code*********************************************/

    就是这样会出错,总是不能正确创建带有坐标的shp矢量文件,那么问题出在哪里?各种尝试,各种调试都不行,最终实在受不了了,改成了如下代码:

/**********************************key_code***********************************************/

        //原始栅格影响空间参考信息直接赋值给矢量文件

      
                
            InputOGRS.SetProjCS("UTM 17 /WGS84");
            InputOGRS.SetWellKnownGeogCS("WGS84");	 	
            InputOGRS.SetUTM(17);	
poLayer = poDS->CreateLayer((const char*)FireArea_SHPpath, InputOGRS, wkbUnknown, NULL);
/*************************************key_code*********************************************/

这样居然好了,我去,那我看一下区别到底在哪里,于是用txt查看***.prj文件,发现了区别:

1.第二种方法生成的.prj文件信息如下:

PROJCS["WGS_1984_UTM_Zone_51N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

2.第一种方法生成的.prj文件信息如下:

PROJCS["unnamed",PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

很明显,第一种方法生成的.prj文件不含参考坐标系,但是检查代码,调试状态查看*WKT和InputOGRS两个参数内容如下:

PROJCS["unnamed",PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]

对比发现第一种方法写入shp的信息在原始信息基础上少了AUTHORITY["EPSG","9001"]我试着直接修改.prj文件,加上AUTHORITY["EPSG","9001"],信息加上只有依然无法在GIS中显示坐标系,可见,根本原因就是缺少了参考坐标,那么针对UTM的原始图像坐标系,如何原封不动的写入shp文件中呢?难道要判断是否是UTM再提取UTMZone投影带号,再用OGR生成后写入?好麻烦的说,于是试了一下如下方法:

/**********************************key_code***********************************************/

        //原始栅格影响空间参考信息直接赋值给矢量文件

        const char *WKT = pInDataset->GetProjectionRef();	

	if (oSRS.SetFromUserInput(WKT ) != OGRERR_NONE)
	    {	       
		    GDALDestroyDriverManager();
	            return 6;//矢量文件写入参考空间信息失败
	    }	
InputOGRS.SetWellKnownGeogCS("WGS84");	
poLayer = poDS->CreateLayer((const char*)FireArea_SHPpath, InputOGRS, wkbUnknown, NULL);
/*************************************key_code*********************************************/

我只是加了一句:

InputOGRS.SetWellKnownGeogCS("WGS84");	

    问题就解决了,为什么呢,因为我的项目都是84的,如果需要其他参考系,应该也可以参考这种方法吧,希望能够对大家有所帮助,不要像我一次次试,一次次失败后才发现根本原因;其实,这也反应了一个问题,那就是分析问题解决问题的思路,以及问题定位的思路;

转载自:https://blog.csdn.net/xujie126/article/details/79557500