[gdal-dev] Question to warping images using python-bindings

[Even Rouault]

Le mercredi 18 février 2015 15:09:14, Jacob Mendt a écrit :
> Dear Gdal-Team and List,
> 
> 
> 
> At first I want to thank you for the great work you guys are doing with
> GDAL. I really enjoy using this tools. I am especially doing rectification
> of unreferenced raster images with the gdal tools. At the beginning I did
> the rectification with the gdal_translate and gdalwarp tools. This works
> pretty well and allowed me to use different order polynomials and a plate
> spline transformer for my image transformation. Also it was possible for me
> to use different resampling algorithm for my images.
> 
> 
> 
> The problem was that I always had to clip the input image (using pixel
> coordinates)  first, save the result in a temporally dataset - first I/O
> process - and doing than the transformation using gdal_translate (add gcps)
> and gdalwarp (transformation and resampling) - second I/O process. In my
> use case the most time the processes spent was for the I/O operations. So
> I decided to use the python-bindings for processing, in hoping to do some
> operations in memory.
> 
> 
> 
> This works great for doing a affine transformation
> (dataset.SetGeotransform) without resampling and gave a big performance
> improvement. The improvement came through reading the file once in memory,
> doing there the clipping and adding of the transformation matrix and then
> writing this to disk. So I got only once to read the image from the disk
> and once to write the clipped image, together with the transformation
> matrix, back to disk. But I couldn't figure out how to do for example a
> plate spline transformation or an order 2 polynomial transformation. After
> further investigation it seems like this is not possible with the python
> bindings, are I am right with this? Also it seems that the python bindings
> doesn't allow to do a resampling of an image without using the
> gdal.ReprojectImage method, what from my understand would make sense in my
> case or?

Jacob,

it seems you're right. The C/C++ GDALReprojectImage() has a final argument 
GDALWarpOptions *psOptions whose char **papszWarpOptions member could be used 
to specify the transformation method with METHOD=TPS or MAX_GCP_ORDER=2. 
Unfortunately the GDALWarpOptions structure hasn't been mapped yet into the 
SWIG bindings, so a NULL pointer is currently used for psOptions (2nd order 
polynomial should be used by default if there are 6 or more GCPs defined and no 
geotransform)
I think it shouldn't be too hard to map the structure (at least the relevant 
parts of it, that aren't overriden by GDALReprojectImage()) into SWIG and make 
it usable by gdal.ReprojectImage().

Even


-- 
Spatialys - Geospatial professional services
http://www.spatialys.com