[gdal-dev] Converting between row/col numbers of two datasets

Knut-Frode Dagestad knutfrodesoppel at hotmail.com
Wed Jan 18 11:01:14 EST 2012


Hi group,

In searching for the reason for the bug reported in ticket #4442 [1], I 
have made some tests with gdal.Transformer in Python:

A transformer defined generally as:
T = gdal.Transformer(src_ds, dst_ds, [''])
can be used to convert row/col from src-image to row/col of dst-image. 
Behind the scene, according to GDAL documentation [2], this consists of 
3 steps:

(source row/col) --> (source coordinates) --> (destination coordinates) 
--> (destination row/col)

My tests show that this works fine if the src-image has GCPs only and 
the dst-image has a geotransform, but not vice versa.
In the latter case, the output is not row/col of the GCP-image, but 
rather the coordinates (lon and lat). I.e., the last step of the chain 
above has not been performed. Enforcing the same projection for the two 
datasets makes no difference.
This should explain why in the tests of ticket #4442 the 
geotransform-image is placed in the upper left corner of the output 
GCP-image: the lon/lat-values are mistaken as destination row/col-values.

However, it is also possible to define a transformer which converts 
(back and forth) between row/column and coordinates of either the src- 
or dst-images:
T_src = gdal.Transformer(src_ds, None, [''])
T_dst = gdal.Transformer(dst_ds, None, [''])
These correspond to the first and last steps of the chain above. The 
interesting thing is that combining these it works well to calculate 
row/col of the GCP-image from row/col of the geotransform-image.

Can anyone explain why then the direct transformation does not work? I 
suspect a small fix could possibly solve this problem?

Best regards from Knut-Frode


[1] http://trac.osgeo.org/gdal/ticket/4442
[2] http://www.gdal.org/gdal__alg_8h.html#7671696d085085a0bfba3c3df9ffcc0a


The script below can be used to reproduce this:

GCPfile = 'MER_RR__2PNPDK20030809_100609_000022512018_00466_07534_3898.N1'
# http://envisat.esa.int/services/sample_products/meris/RR/L2/

geotransformfile = 'europe.tif'# 
http://www.eea.europa.eu/data-and-maps/data/world-digital-elevation-model-etopo5/zipped-dem-geotiff-raster-geographic-tag-image-file
-format-raster-data/

# ----- Open source and destination datasets -----#
if 0:
         # This works
         src_ds = gdal.Open(GCPfile)
         dst_ds = gdal.Open(geotransformfile)
else:
         # This does not work
         src_ds = gdal.Open(geotransformfile)
         dst_ds = gdal.Open(GCPfile)

# ----- Make transformers -----#

# To tranform from pixel/line of src-image to pixel/line on dst-image:
T_full = gdal.Transformer(src_ds, dst_ds, [''])

# To tranform between pixel/line and coordinates on src-image
T_src = gdal.Transformer(src_ds, None, [''])

# To tranform between pixel/line and coordinates on dst-image
T_dst = gdal.Transformer(dst_ds, None, [''])

# ----- Test transformers -----#
src_rowcol = (100, 100)
coords = T_src.TransformPoint(0, src_rowcol[0], src_rowcol[1])[1]
print 'Pixel ' + str(src_rowcol) + ' in src image has coordinates: ' + 
str(coords)

dst_rowcol = T_dst.TransformPoint(1, coords[0], coords[1])[1];
print 'These lon/lat-values corresponds to pixels in destination image: 
' + str(dst_rowcol)

print 'Direct transformation: ' + str(T_full.TransformPoint(0, 
src_rowcol[0], src_rowcol[1]))





More information about the gdal-dev mailing list