Using the Geospatial Data Abstraction Library (GDAL) in Ubuntu

For a project, I need to process SRTM elevation data. GDAL is a software library that makes it easy to read and write raster geospatial data formats. I decided to give it a try since SRTM data uses the ESRI grid format supported by the library.

1. Install the the GDAL libraries and utilities

$ apt-get install libgdal1-1.5.0 gdal-bin

2. Sample code. The SRTM data file path is in the variable pszFilename. The code is based on the tutorial.

#include <stdio.h>
#include <gdal.h>>
int main()
{
    GDALDatasetH  hDataset;
    char *pszFilename = "/media/DATA/srtm_61_10.asc";
    GDALDriverH   hDriver;
    double        adfGeoTransform[6];
    GDALRasterBandH hBand;
    int             nBlockXSize, nBlockYSize;
    int             bGotMin, bGotMax;
    double          adfMinMax[2];
    float *pafScanline;
    int   nXSize;
    int i;

    GDALAllRegister();
    hDataset = GDALOpen( pszFilename, GA_ReadOnly );
    if( hDataset == NULL )
    {
        printf("Cannot open file.");
    }
    hDriver = GDALGetDatasetDriver( hDataset );
    printf( "Driver: %s/%s\n",
            GDALGetDriverShortName( hDriver ),
            GDALGetDriverLongName( hDriver ) );
    printf( "Size is %dx%dx%d\n",
            GDALGetRasterXSize( hDataset ), 
            GDALGetRasterYSize( hDataset ),
            GDALGetRasterCount( hDataset ) );
    if( GDALGetProjectionRef( hDataset ) != NULL )
        printf( "Projection is `%s'\n", GDALGetProjectionRef( hDataset ) );
    if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
    {
        printf( "Origin = (%.6f,%.6f)\n",
                adfGeoTransform[0], adfGeoTransform[3] );
        printf( "Pixel Size = (%.6f,%.6f)\n",
                adfGeoTransform[1], adfGeoTransform[5] );
    }
    hBand = GDALGetRasterBand( hDataset, 1 );
    GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
    printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
            nBlockXSize, nBlockYSize,
            GDALGetDataTypeName(GDALGetRasterDataType(hBand)),
            GDALGetColorInterpretationName(
            GDALGetRasterColorInterpretation(hBand)) 
    );
    adfMinMax[0] = GDALGetRasterMinimum( hBand, &bGotMin );
    adfMinMax[1] = GDALGetRasterMaximum( hBand, &bGotMax );
    if( ! (bGotMin && bGotMax) )
      GDALComputeRasterMinMax( hBand, TRUE, adfMinMax );
    printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
    if( GDALGetOverviewCount(hBand) > 0 )
      printf( "Band has %d overviews.\n", GDALGetOverviewCount(hBand));
    if( GDALGetRasterColorTable( hBand ) != NULL )
      printf( "Band has a color table with %d entries.\n", 
                 GDALGetColorEntryCount(
                 GDALGetRasterColorTable( hBand ) ) 
      );
    nXSize = GDALGetRasterBandXSize( hBand );
    pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize);
    GDALRasterIO( hBand, GF_Read, 0, 0, nXSize, 1, 
                    pafScanline, nXSize, 1, GDT_Float32, 
                     0, 0 );
    for (i=0;i<20;i++){
      printf("%f\n", pafScanline[i]);
    }
}

3. Build the example

$ gcc -lgdal1.5.0 -I/usr/include/gdal -o testgdal.exe testgdal.c

4. Run

$./testgdal.exe
Driver: AAIGrid/Arc/Info ASCII Grid
Size is 6001x6001x1
Projection is `GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]]'
Origin = (119.999583,15.000417)
Pixel Size = (0.000833,-0.000833)
Block=6001x1 Type=Int16, ColorInterp=Undefined
Min=-27.000d, Max=2375.000