Search code examples
c#rastergdalprojection

GDAL ReadRaster get pixel block map extent

GDAL's user guide on geotransform between pixel row/col to map coordindates says to use the following parameter array:

GT[0] x-coordinate of the upper-left corner of the upper-left pixel.
GT[1] w-e pixel resolution / pixel width.
GT[2] row rotation (typically zero).
GT[3] y-coordinate of the upper-left corner of the upper-left pixel.
GT[4] column rotation (typically zero).
GT[5] n-s pixel resolution / pixel height (negative value for a north-up image)

So, if I have a pixel block that is 10 columns to the right and 20 rows down from the upper left corner's coordinate, i.e. GT[0], GT[3], how to calculate the extent of the 10 by 20 block of pixels?

I would think the pixel block's extent would consist of 4-corner's map (x, y), for example:

the upper left vertex (x, y) would be:
  (x_map = GT[0], 
   y_map = GT[3])

the upper right vertex (x, y) would be:
  (x_map = GT[0] + 10 * GT[1] + 20 * GT[2], 
   y_map = GT[3])

the lower right vertex (x, y) would be:
  (x_map = GT[0] + 10 * GT[1] + 20 * GT[2],
   y_map = GT[3] + 10 * GT[4] + 20 * GT[5])

the lower left vertex (x, y) would be:
  (x_map = GT[0],
   y_map = GT[3] + 10 * GT[4] + 20 * GT[5])

here is a real example: enter image description here does this look alright? is there a GDAL way of getting these directly? what if the grid is really big, meaning the original upper left corner could be really far away from the pixel block I am trying to get the extent from, wouldn't the error accumulate in calculated vertex's map coordinate?

Solution

  • You can use the ApplyGeoTransform function to have GDAL do this. I'm using the Python API, but it should be the same/similar in C#:

    from osgeo import gdal

gt = (-124.5, 9.26e-5, 0, 50.0017, 0, -9.26e-5)

# upper left
ulx, uly = gdal.ApplyGeoTransform(gt, 0, 0)

# lower right
lrx, lry = gdal.ApplyGeoTransform(gt, 10+1, 20+1)

    Which gives (-124.5, 50.0017, -124.4989814, 49.9997554) as the extent (for ulx, uly, lrx, lry). With this tiny resolution there might be a bit of rounding happening in the values above. Depending on how you define your block, you might need to add 1 pixel for the lower right.

    Note that it gives the upper left location of a pixel, not the center. To get the center you can add half a pixel (eg gdal.ApplyGeoTransform(gt, 10+0.5, 20+0.5). That's different to how for example formats like NetCDF often store coordinates.

    There's is an accompanying function that inverts the geotransform, which allows you to calculate the reverse (map to pixel).

    gt_inv = gdal.InvGeoTransform(gt)
ulx_px, uly_px = gdal.ApplyGeoTransform(gt_inv, lrx, lry)

    Using the same example gives me (10.00000000023283, 20.0) (for ulx_px, uly_px). The function returns fractional pixel coordinates, so you might need to cast them to integer depending on the application.

    If your resolution in the geotransform is not exact, it will indeed accumulate an error. That's inevitable with finite precision. So choosing the extent and resolution a bit careful if possible can help. Storing a known extent in the metadata could be a workaround if you really need a "clean" extent without having to rely on the geotransform, but that would be unconventional.