5

This question already has an answer here:

I'm now stuck with writing a piece of code to look up individual pixel values on a raster image.

This piece of code mainly executes two tasks: (1) converting world coordinates to raster coordinates (I use a package named "Affine" to do that); (2) then computing pixel offsets and figuring out which column and row the pixel is located at.

The major problem I'm facing now is that the code didn't return the correct column and row that the pixel is located at.

import sys, gdal, numpy
from affine import Affine

# coordinates to get pixel values for
lons = [-122.265373, -122.429139, -68.123095]
lats = [37.873090, 37.783640, 44.981935]

# open the raster file
ds = gdal.Open('data/nlcd_compressed.tif')
if ds is None:
    print 'Could not open the raster file'
    sys.exit(1)

# get raster size
cols = ds.RasterXSize #161190 
rows = ds.RasterYSize #104424
bands = ds.RasterCount #1

# get georeference info
transform = ds.GetGeoTransform() # (-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0)
xOrigin = transform[0] # -2493045.0
yOrigin = transform[3] # 3310005.0
pixelWidth = transform[1] # 30.0
pixelHeight = transform[5] # -30.0

# convert global coords to raster coords
aff = Affine.from_gdal(-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0)
x_coords, y_coords = aff*(numpy.array(lons), numpy.array(lats))
# converted coordinates
xValues = [-2496712.96119, -2496717.87417, -2495088.69285]
yValues = [3308868.8073, 3308871.4908, 3308655.54195]


band = ds.GetRasterBand(1) # 1-based index
#bandtype = gdal.GetDataTypeName(band.DataType) #Byte
#data = band.ReadAsArray(0, 0, cols, rows)
#value = data [122, 37] #127

data = band.ReadAsArray(xOffset)

# loop through the coordinates
for i in range(len(lons)):
    x = xValues[i]
    y = yValues[i]
    # compute pixel offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)
    # get individual pixel values
    data = band.ReadAsArray(xOffset, yOffset, 1, 1)
    value = data[0, 0]
    print value

marked as duplicate by user2856, John Powell, Andre Silva, nmtoken, whyzar Jan 22 '17 at 17:29

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • Have you checked gdal2xyz which converters your entire raster into lat, long, value and gdallocationinfo which will tell the raster value at any latlong or pixel xy. – GeoSpatialEarth.in Jan 22 '17 at 5:47
  • I rewrote your code again to avoid unnecessary lines and to express coordinates as points. It works perfectly. – xunilk Jan 22 '17 at 10:02
5

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;
}
4

I rewrote your code again to avoid unnecessary lines and to express coordinates as points. I tried out it with my particular raster (utah_demUTM2.tif) and it works perfectly (see next image).

from osgeo import gdal
import sys

# coordinates to get pixel values for (as tuples of points)
points = [(401229.289973, 4466415.48331)]

# open the raster file
ds = gdal.Open('/home/zeito/pyqgis_data/utah_demUTM2.tif')

if ds is None:
    print 'Could not open the raster file'
    sys.exit(1)
else:
    print 'The raster file was opened satisfactorily'

# get georeference info
transform = ds.GetGeoTransform() # (-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0)
xOrigin = transform[0] # -2493045.0
yOrigin = transform[3] # 3310005.0
pixelWidth = transform[1] # 30.0
pixelHeight = transform[5] # -30.0

band = ds.GetRasterBand(1) # 1-based index

data = band.ReadAsArray()

# loop through the coordinates
for point in points:
    x = point[0]
    y = point[1]

    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)
    print xOffset
    print yOffset
    # get individual pixel values
    value = data[yOffset][xOffset]
    print value

After running the code at Python Console of QGIS, I corroborated printed results (value, xOffset, yOffset) with raster value, by using Value Tool plugin, and they matched.

enter image description here

2

Alternatively you can use rasterio, which has a sample() method.

import rasterio

lons = [-122.265373, -122.429139, -68.123095]
lats = [37.873090, 37.783640, 44.981935]

with rasterio.open('/path/to/file.tif') as src:
    for val in src.sample(zip(lons, lats)):
        print val

Use a list comprehension if you want a list of values

with rasterio.open('/path/to/file.tif') as src:
    vals = [x for x in src.sample(zip(lons, lats))]

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