5

I am trying to gather information from a *.tif file using the GDAL libraries. When I open the *.tif file using both the Python API, and C++ API, I encouter two completely different answers. Here are the two blocks of code, and their terminal outputs:

Python

import gdal
from gdalconst import GA_ReadOnly

def main():
    data_set = gdal.Open("file_name.tif", GA_ReadOnly)
    print data_set.GetProjectionRef()
    print 'Driver: ', data_set.GetDriver().LongName

if __name__ == '__main__':
    main()

Output

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

Driver: GeoTIFF

C++

#include "gdal_priv.h"

int main() {
  GDALDataset  *poDataset;

  GDALAllRegister();

  poDataset = (GDALDataset *) GDALOpen("file_name.tif", GA_ReadOnly);

  if(poDataset == NULL)
    return 0;

  if(poDataset->GetProjectionRef() != NULL)
    printf("Projection is `%s'\n", poDataset->GetProjectionRef());

  printf(
    "Driver: %s/%s\n", poDataset->GetDriver()->GetDescription(), 
    poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));

  return 0;
}

Output

Projection is

`PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32613"]]'
Driver: GTiff/GeoTIFF

If it helps in the answer. I used the terminal to install the C libraries, and I am using Anaconda to manage my Python libraries.

  • and the result of gdalinfo file_name.tif ? – gene Apr 7 '16 at 16:17
  • @gene: Driver: GTiff/GeoTIFF Files: file_name.tif Size is 9500, 8541 Coordinate System is: PROJCS["unnamed", PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-105], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]]] Origin = (479274.000000000000000,4426986.000000000000000) Pixel Size = (1.00000000000000,-1.00000000000000) Image Structure Metadata: INTERLEAVE=BAND Band 1 Block=9500x1 Type=Float32, ColorInterp=Gray – Dylan Hicks Apr 7 '16 at 16:25
  • I don't think the 2 different outputs are completely different. They seem to be functionally equivalent. – Jesse Crocker Apr 7 '16 at 16:41
  • @JesseCrocker You are right. I should rephrase. The output of the Python code seems to be a subset of the output of the C++ code. – Dylan Hicks Apr 7 '16 at 16:44
  • Make sure GDAL_DATA is properly set for both instances. – user10353 Apr 11 '16 at 22:26
1

If the output from gdalinfo file_name.tif is to be believed, it looks like the C++ implementation assumed your datum was WGS84 and filled out the rest of the WKT for you, although there shouldn't be an entry for both PROJCS and GEOGCS.

It's hard to tell for certain without digging into the gdal source code and libraries.

1

You must understand the logic of the Python version of GDAL. With a raster with the same projection (EPSG:3857)

from osgeo import gdal
data_set = gdal.Open("test.tif")
# Getting the Spatial Reference System (returned in Well Known Text(WKT) format)
projInfo = data_set .GetProjection()
projInfo
'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],AUTHORITY["EPSG","3857"]]'
data_set.GetDriver().LongName
'GeoTIFF'

We can also use the osr module to create a spatial reference object

from osgeo import osr
Ref = osr.SpatialReference()
Ref.ImportFromWkt(projInfo)
# Proj.4 format
RefProj = Ref.ExportToProj4()
print Ref
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
AUTHORITY["EPSG","3857"]]
print RefProj
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.