[GIS] How to get coordinate system info from shapefile using GDAL/OGR

ccoordinate systemgdalogrshapefile

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?

Best Answer

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())
Related Question