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
3. Build the example
1. Install the the GDAL libraries and utilities
$ apt-get install libgdal1-1.5.0 gdal-bin2. 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.c4. 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