[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