[C/C++] GIS 따라해보기 – gdal/GEOS SHP 파일 조작

우선 shp 파일을 조작 해 보기 위해서는 해당 파일을 구해야 하는 데 구글링 해보니 다음과 같은 사이트가 있었습니다. 해당 사이트에서 파일을 다운로드 받아봐도 될 것 같습니다

 

공짜로 데이터얻는 데 광고도 함 넣어주고… 데이터 받아서 QGIS에서 열어보니, 강인듯…보이네요

 

미국 해양 대기청사이트에서도 데이터를 아래와 구할 수 있습니다.

 

NOAA 라….

 

튜토리얼에 따라서위의 폴리라인 SHP을 읽어 들이는 코드를 작성 한다면 다음 사이트를 보고 작성 해도 될 것입니다.

https://gdal.org/tutorials/vector_api_tut.html

 

Vector API tutorial — GDAL documentation

 

#include "ogrsf_frmts.h"

int main()
{
	GDALAllRegister();
	GDALDataset *poDS;

	poDS = (GDALDataset*) GDALOpenEx( "KOR_wat/KOR_water_lines_dcw.shp", GDAL_OF_VECTOR, NULL, NULL, NULL );

	if( poDS == NULL )
	{
	   printf( "Open failed.\n" );
	   exit( 1 );
	}

	OGRLayer *poLayer;
	poLayer = poDS->GetLayerByName( "KOR_water_lines_dcw" );

	OGRFeature *poFeature;
	poLayer->ResetReading();

	while( (poFeature = poLayer->GetNextFeature()) != NULL )
	{

		OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
		int iField;

		for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
		{

			OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
			
			if( poFieldDefn->GetType() == OFTInteger )
				printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
			else if( poFieldDefn->GetType() == OFTInteger64 )
				printf( CPL_FRMT_GIB ",", poFeature->GetFieldAsInteger64( iField ) );
			else if( poFieldDefn->GetType() == OFTReal )
				printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
			else if( poFieldDefn->GetType() == OFTString )
				printf( "%s,", poFeature->GetFieldAsString(iField) );
			else
				printf( "%s,", poFeature->GetFieldAsString(iField) );
				//printf("\r\n");
		}

		char *pszWKT = NULL;

		OGRGeometry *poGeometry, *poGeometry2;

		poGeometry = poFeature->GetGeometryRef();
		poGeometry->exportToWkt(&pszWKT);

		printf( "%s\n", pszWKT );

		CPLFree(pszWKT);
		OGRFeature::DestroyFeature( poFeature );

	}

	GDALClose( poDS );

}

 

 

우선 SHP 파일을 열기 전에 해야 하는 작업은 해당 SHP 파일을 읽어들이기 위한

드라이버를 생성하는 것인데, 여기서는 가용한 모든 드라이버를 열었다고 보면 됩니다.

GDALAllRegister();

 

그런 다음 GDALOpenEx 함수를 이용해서 SHP 파일을 연다음, GDALDataset 클래스에 저장 해주고,

GDALDataset *poDS;
poDS 
=(GDALDataset*) GDALOpenEx( "KOR_wat/KOR_water_lines_dcw.shp",GDAL_OF_VECTOR, NULL, NULL, NULL );

 

해당 Feature의 속성정보를 얻기 위해서는 OGRFeatureDefn 클래스의참조 값을 얻어서 다음과

같이 코딩하여주면 되는 듯 합니다.

 

다음 코드로 dbf 속성값을 볼 수 있습니다:

OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();

for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
{

	OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
	
	if( poFieldDefn->GetType() == OFTInteger )
		printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
	else if( poFieldDefn->GetType() == OFTInteger64 )
		printf( CPL_FRMT_GIB ",", poFeature->GetFieldAsInteger64( iField ) );
	else if( poFieldDefn->GetType() == OFTReal )
		printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
	else if( poFieldDefn->GetType() == OFTString )
		printf( "%s,", poFeature->GetFieldAsString(iField) );
	else
		printf( "%s,", poFeature->GetFieldAsString(iField) );
		//printf("\r\n");
}

 

그냥 뿌리면아래와 같이 떡이 되는 듯 한데 포맷이 필요 할 듯 하네요...

 

그 이후에 공간자료의 좌표속성을 WKT 즉 Well Known Text 형태로뿌려 보았다.

공간자료 좌표를 뿌리는 것은 OGRGeometry 의 exportToWkt 함수를 사용하였으며, 다음과 같이 사용하면 됩니다.

poGeometry = poFeature->GetGeometryRef();
poGeometry->exportToWkt(&pszWKT);
printf( "%s\n", pszWKT );
CPLFree(pszWKT);

 

 

문서 상에서도 말했듯이, 생성된 WKT 문자열 해제는 꼭 해줘야 하며, 다음과 같이 CPLFree 함수를 사용한다고 합니다.

CPLFree(pszWKT);

 

출력화면:

여타 폴리곤을열어서 속성을 확인 하는 작업도 마찬가지이며, 예제는 포인트 자료를 열어 보는 것으로 진행 하게 됩니다.

 

여기서 소스코드를 gcc로 컴파일 하고 MinGW 쉘 상에서 구동하면 안되고, 따로 윈도우 용 cmd 창을 띄워서 실행 하도록 하였는 데, 이렇게 하면 아래와 같이 경로 설정이 필요 합니다.

SET PATH=C:\DEV\COMP\msys32\home\tobee\bin\gdal\bin;
C:\DEV\COMP\msys32\home\tobee\bin\geos\bin;
C:\DEV\COMP\msys32\mingw32\bin;
C:\DEV\COMP\msys32\usr\bin

 

경로 설정은 상황에 따라 틀린데, 유의 할 점이라면 GEOS/GDAL/MinGW32 의 bin 디렉토리를 잡아주고, mingw 용 libexpat-1.dll이 있음에도 불구하고, msys용 expat 즉 msys-libexpat-1.dll을 찾고 있네요. 뭔가 빌드가 잘 못 된듯 한데 넘어가기로 합니다...

 

이상으로 LineString을 구성된 SHP 파일을 열어보고 해당 SHP 파일의 속성인 dbf 파일 보기와 WKT 형태의 공간자료를 뽑아 보았습니다.

 

이상.