GDAL has all tools on board to read the data. The steps howto read data can be found in http://www.gdal.org/gdal_tutorial.html section read data. The routines GDALRasterIO or in python band.ReadRaster at least will do the job in conjunction with the inverse-pixel-to world-transformation given by GDALGetGeoTransform (parameter of the transform) and calcWorldToPixel.
Here is a step by step procedure (unfortunally only for the C-API) for an orientation, I've written last year. In my test setup I call:
./read-data test.tif 7.8146154 54.2904329
TRANSFORM:
X = 7.812837e+00 + 2.379156e-07 * COL + 0.000000e+00 * ROW
Y = 5.429176e+01 + 0.000000e+00 * COL + -2.379156e-07 * ROW
LAYOUT: WxHxB is 14946x11192x3 Bands found!
FILE: test.tif WORLD: 7.814615e+00 5.429043e+01 PIXEL: 7473 5596
BAND: 1 TYPE: UInt16 KEY: 2 SIZE: 16
BAND: 2 TYPE: UInt16 KEY: 2 SIZE: 16
BAND: 3 TYPE: UInt16 KEY: 2 SIZE: 16
The result for an 16 Bit RGB is:
14946x11192=[ 8906 18372 23424 ]
Here is the documented code:
// -------------------------------------------------------------
// Example code how to fetch data from an GDAL raster source
// -------------------------------------------------------------
// gcc -lgdal -lm -I /usr/include/gdal -std=c99 read-pixel.c -o read-pixel
// -------------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdint.h>
#include "gdal.h"
#include "cpl_conv.h"
// Maschine ZERO
#define DBL_EPSILON 2.2204460492503131e-16
#define LINE_WISE 0
#define SBUF_SIZE 16
#define BBUF_SIZE 96
I use the both transformation forms based on the parameter arrays given by GDALGetGeoTransform and calcWorldToPixel is the interesting one.
// --------------------------------------------------------------
// Transformation World to Pixel
// trfm - Parameter array read from gdal raster dataset
// x, y - World postion
// col, row - Pixel postions (return values)
// return true if the calculation is valid and 0 if there is a
// division by zero
// -------------------------------------------------------------
int calcWorldToPixel(double *trfm,
double x, double y,
long *col , long *row) {
double div = (trfm[2]*trfm[4]-trfm[1]*trfm[5]);
if (div<DBL_EPSILON*2) return 0;
double dcol = -(trfm[2]*(trfm[3]-y)+trfm[5]*x-trfm[0]*trfm[5])/div;
double drow = (trfm[1]*(trfm[3]-y)+trfm[4]*x-trfm[0]*trfm[4])/div;
*col = round(dcol); *row = round(drow);
return 1;
}
// --------------------------------------------------------------
// Transformation pixel to world
// trfm - Parameter array read from gdal raster dataset
// col, row - Pixel postions
// x, y - World postion (return values)
// -------------------------------------------------------------
void calcPixelToWorld(double *trfm,
double col, double row,
double *x , double *y) {
*x = trfm[0] + trfm[1] * col + trfm[2] * row;
*y = trfm[3] + trfm[4] * col + trfm[5] * row;
}
The main program looks like this:
// -------------------------------------------------------------
// Main with parameter:
// arg[1] File name of the raster source
// arg[2] Longitude position
// arg[3] Latitude position
// success if return is zero ..surprise
// -------------------------------------------------------------
int main(int argc, char* argv[]) {
// Test number of arguments
if (argc<4) {
printf("Usage %s pix file lon lat\n",argv[0]);
return 1;
}
// Register GDAL drivers
GDALAllRegister();
// Declare the dataset handler
GDALDatasetH hDataset;
// and the affine transformation
double trfm[6];
// Geotiff file argument 1
char *pFileName = argv[1];
// Read LON, LAT from commandline args
double lon = 0.0; double lat = 0.0;
if (! sscanf(argv[2],"%lf",&lon) ) {
printf("Invalid numerical lon value for %s\n!",argv[2]);
return 10;
};
if (! sscanf(argv[3],"%lf",&lat) ) {
printf("Invalid numerical lat value for %s\n!",argv[3]);
return 20;
};
// Declare the pixel vars
long col = -1; long row = -1;
OK here its starts with the investigation of contained data to calculate the stuff:
// Try to open the raster data set
hDataset = GDALOpen( pFileName, GA_ReadOnly );
if( hDataset == NULL ) {
printf("Cannot open file %s !\n", pFileName);
return 30;
}
// Read transform from raster
if( GDALGetGeoTransform( hDataset, trfm ) == CE_None ) {
printf("TRANSFORM: \n");
printf(" X = %e + %e * COL + %e * ROW\n", trfm[0], trfm[1], trfm[2] );
printf(" Y = %e + %e * COL + %e * ROW\n", trfm[3], trfm[4], trfm[5] );
} else {
printf("Missing transformation in %s !\n",pFileName);
GDALClose(hDataset);
return 40;
}
Get the image layout and the data type:
int imgWidth = GDALGetRasterXSize( hDataset );
int imgHeight = GDALGetRasterYSize( hDataset );
int numBands = GDALGetRasterCount (hDataset );
if ( numBands <1 ) {
printf("Missing band info in %s",pFileName);
GDALClose(hDataset);
return 50;
} else {
printf("LAYOUT WxHxB is %dx%dx%d Bands found!\n",
imgWidth, imgHeight, numBands);
}
// Determin image dimension and type
GDALRasterBandH hBand[numBands];
GDALDataType hType[numBands];
int hSize[numBands];
// Determin band type
for(int b=0 ; b<numBands; b++) {
hBand[b] = GDALGetRasterBand( hDataset, b+1 );
hType[b] = GDALGetRasterDataType(hBand[b]);
hSize[b] = GDALGetDataTypeSize(hType[b]);
printf("BAND: %d TYPE: %s KEY: %d SIZE: %d\n", b+1,
GDALGetDataTypeName(hType[b]), hType[b], hSize[b]);
}
Now you could calculate the pixel pos using the calcWorldToPixel(...):
// Calc pixel position
int res = calcWorldToPixel(trfm,lon,lat,&col,&row);
if ( !res || row<0 || col <0 ||
col >= imgWidth-1 || row >= imgHeight-1) {
printf("Cannot calculate data for FILE: %s POS: %e %e!\n",
pFileName, lon ,lat);
GDALClose(hDataset);
return 60;
} else {
printf("FILE: %s WORLD: %e %e PIXEL: %ld %ld \n",
pFileName, lon ,lat, col, row);
}
And retrieve the data with GDALRasterIO in Python (band.ReadRaster(..)). In C I've to distinguish the data type of the bands only done here for GDT_UInt16.
// Collect the data sets
for(int b=0 ; b<numBands; b++) {
// Handle the data type
switch(hType[b]) {
// More types here
// case GDT_Byte: {
// }
case GDT_UInt16: {
Work with a line buffer
if (LINE_WISE) {
// Allocate a line buffer
uint16_t *pData = (uint16_t *) CPLMalloc(hSize[b]*imgWidth);
// Read it
GDALRasterIO( hBand[b],
GF_Read, // Read from band
0, row, // Offs X, Offs Y
imgWidth, 1, // Exact one line
pData, // The target buffer
imgWidth, 1, // Size of the target buffer
hType[b], // Buffer type
0, 0 ); // Array strides X/Y dir
// Format and write it to the collector
sprintf(buffer," %d ", pData[col]);
strncat(collect, buffer, BBUF_SIZE);
// Release it
CPLFree(pData);
} else {
} else {
Or a single cell wich is slow for frequent requests to the file
// Declare a single cell
uint16_t pData;
// Read it
GDALRasterIO( hBand[b],
GF_Read, // Read from band
col, row, // Offs X, Offs Y
1, 1, // Exact one value
&pData, // The target buffer
1, 1, // Size of the target buffer
hType[b], // Buffer type
0, 0 ); // // Array strides X/Y dir
// Collect it
sprintf(buffer," %d ",pData);
strncat(collect, buffer, BBUF_SIZE);
}
break;
}
default: {
printf("Band type is not suported\n");
return 70;
break;
}
}
}
At least print the result.
printf("%s] \n", collect);
GDALClose(hDataset);
return 0;
}
