How to get coordinate system info from shapefile using GDAL/OGR?

Question

GDAL version 2.1.1

If I open a geospatial raster file, such as a GeoTIFF or DTED file, with GDAL, then I can get the coordinate system info for that file, as a string in OpenGIS WKT format. For example:

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]

The C++ code I use to get that info is simple and straightforward.

string path = "somewhere.dt1";
unsigned int openFlags = GA_ReadOnly;
GDALDataset * pDataset = static_cast<GDALDataset *>( GDALOpenEx(
    path.c_str(),
    openFlags,
    NULL, NULL, NULL ) );

string projWkt = pDataset->GetProjectionRef();

OGRSpatialReference srs( projWkt.c_str() );
srs.dumpReadable();

However, when I open a geospatial vector file, such as a shapefile, with GDAL, I get only a blank string for the coordinate system info. The only difference in the C++ code I use is that I must also provide the GDAL_OF_VECTOR flag to GDALOpenEx().

string path = "elsewhere.shp";
unsigned int openFlags = GA_ReadOnly | GDAL_OF_VECTOR;
GDALDataset * pDataset = static_cast<GDALDataset *>( GDALOpenEx(
    path.c_str(),
    openFlags,
    NULL, NULL, NULL ) );

The shapefile I am opening does have a valid coordinate system. I can confirm that by opening the shapefile in ArcMap or ArcCatalog. Both of those apps are able to tell me the coordinate system and projection info.

So, why am I not able to fetch it via this GDAL API?

Is there a different mechanism for vector files?

Answer 1

Is there a different mechanism for vector files?

Yes. Although the vector and raster drivers were merged in GDAL 2x, you still need to know which methods apply to rasters and which are for vectors. Confusingly, vector datasources now expose methods that only apply to rasters...

The GetProjectionRef method you are using applies to raster datasources and you need to use the vector method GetSpatialRef (which returns a spatial reference object, not WKT...).

This works for vector (python code):

from osgeo import gdal, ogr

f = r'/path/to/shapefile.shp'

ds = gdal.OpenEx(f, gdal.OF_VECTOR) #gdal.OF_VECTOR is not required, but speeds things up as only vector drivers are tried
print('"%s"'%ds.GetLayer().GetSpatialRef().ExportToWkt())

Answer 2

Based on the answer from @Luke, the following is my C++ implementation.

I detect whether the file is raster or vector, using a utility function mentioned in another post also answered by @Luke.

string path = "elsewhere.shp";
Format::Enum format = GetFileFormat( path );

I open the file.

GDALDataset * pDataset = NULL;
if ( format == Format::RASTER )
{
  unsigned int openFlags = GA_ReadOnly;
  pDataset = static_cast<GDALDataset *>( GDALOpenEx( inPath.c_str(), openFlags, NULL, NULL, NULL ) );
}
else if ( format == Format::VECTOR )
{
  unsigned int openFlags = GA_ReadOnly | GDAL_OF_VECTOR;
  pDataset = static_cast<GDALDataset *>( GDALOpenEx( inPath.c_str(), openFlags, NULL, NULL, NULL ) );
}

I get the spatial reference system for the file.

OGRSpatialReference srs;
if ( format == Format::RASTER )
{
  char const * szWkt = pDataset->GetProjectionRef();
  srs.importFromWkt( (char **)( &szWkt ) );
}
else if ( format == Format::VECTOR )
{
  // Assume just one vector layer in the dataset for now.
  OGRSpatialReference * pOrigSrs = pDataset->GetLayer( 0 )->GetSpatialRef();
  if ( pOrigSrs )
  {
    srs = *pOrigSrs;
  }
}

Now I can get what I want from the spatial reference system.

if ( srs.IsProjected() )
{
  ...
}