[gdal-dev] Warping onto an image with only GCPs with Python

Fri Jan 13 04:33:59 EST 2012

Dear Brian,

thank you very much for your reply!

The '-geoloc' switch does produce the error about points failing to 
transform but the principal limitation is that I cannot use the '-te' 
switch because I cannot know the extent in advance.

My task is to overlay one image (either projected like e.g. Europe map 
or unprojected like e.g. L1B MERIS or MODIS image) onto another 
UNprojected image (MODIS or MERIS image with GCPs or geolocation 
arrays). And GDAL cannot do that correctly neither with gdalwarp (no 
matter which options I use) nor with Python bindings. It produces an 
output where the source image is place in the top left corner of the 
destination image.

Obviously this problem arises from a bug in the GDAL code and I will try 
to issue a descriptive bug report.

All the best!
Anton

 > Anton
 >
 > EOS_SWATH has 2 geolocation arrays, you might try the -geoloc switch to
 > gdalwarp, otherwise it creates a handfull of gcps from the arrays, if
 > you get a error about points failing to transform, set the output extent
 > with the -te switch.
 >
 > also you could try the swath2grid application in
 > https://lpdaac.usgs.gov/tools/modis_reprojection_tool_swath
 >
 > Brian
 >
 >
 > On Thu, 2012-01-12 at 09:23 +0100, Anton Korosov wrote:
 >> Hello everyone!
 >>
 >> Some time ago I've asked about problems with gdalwarp  (the e-mails
 >> below). I tried to utilize Python binding instead but got the same
 >> results. This brief script overlays the Europe map (GeoTIFF with
 >> geotransform) onto a MERIS images (PDS file with GCPs). But the result
 >> looks wrong again - the map is very small and apperas in the upper left
 >> corner (attached file).
 >>
 >> GCPfile='MER_RR__2PNPDK20030809_100609_000022512018_00466_07534_3898.N1'
 >> geotransformfile = 'europe.tif'
 >> outfile = 'colocated.tif'
 >>
 >> # open input datasets
 >> src_ds = gdal.Open(geotransformfile)
 >> dst_ds = gdal.Open(GCPfile)
 >>
 >> # create writable copy of the destination image
 >> driver = gdal.GetDriverByName( "GTiff" )
 >> out_ds = driver.CreateCopy(outfile, dst_ds, 0 )
 >>
 >> # reproject onto destination image with GCPs
 >> gdal.ReprojectImage(src_ds, out_ds)
 >>
 >> Can you please confirm that this is a bug in GDAL? If so I will issue a
 >> bug report.
 >>
 >> Thank you in advance!
 >> Anton
 >>
 >> PS. If the GCPs in the destination image are replaced with a
 >> GeoTransfrom generated using the gdal.GCPsToGeoTransform() function, the
 >> result looks much better but the image is still rather distorted. Which
 >> is understandable because curvature of a MERIS images (given by GCPs)
 >> cannot be fully explained by few numbers in the geotransform.
 >>
 >>
 >>
 >>
 >> On 12/16/2011 03:56 PM, Anton Korosov wrote:
 >>> Hi!
 >>>
 >>> I'm not sure but probably my question is related to the one below.
 >>> I'm trying to mosaic one MODIS images onto another with GDAL:
 >>>
 >>> gdal_translate
 >>> 
HDF4_EOS:EOS_SWATH:"/Data/sat/GDAL_test/MYD021KM.A2011228.0925.005.2011229003113.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB
 >>> -b 1 /data/modis1.tif
 >>>
 >>> gdal_translate
 >>> 
HDF4_EOS:EOS_SWATH:"/Data/sat/GDAL_test/MOD021KM.A2011229.1125.005.2011229195243.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB
 >>> -b 1 /data/modis2.tif
 >>>
 >>> gdalwarp modis1.tif modis2.tif
 >>>
 >>> Bu the results looks weird - the first image is on top of the 
second one
 >>> but it is very small, distorted and in the top-left corner (though they
 >>> cover approximately the same area).
 >>>
 >>> Do I understand right that such operation would work only when the
 >>> target image georeference is defined by GeoTransform and not by GCPs?
 >>>
 >>> Can I also ask if it is a big work to include ability of GDAL to warp
 >>> onto grid with GCPs (or RPC or geolocation array) ? And if it is 
planned
 >>> in some future?
 >>>
 >>> Thank you!
 >>> Anton
 >>>
 >>> On 12/15/2011 03:47 PM, Knut-Frode Dagestad wrote:
 >>>> Hi,
 >>>>
 >>>> Hróbjartur Þorsteinsson explained to me a near-hidden feature of GDAL,
 >>>> namely to reproject one image onto the coverage of another.
 >>>>
 >>>> Given two images:
 >>>> imageA.tiff
 >>>> imageB.tiff
 >>>>
 >>>> Then an empty image can be created from image A with e.g.
 >>>> $ gdal_translate -ot Float32 -scale 0 0 999 999 -a_nodata 999
 >>>> imageA.tiff domainA.tiff
 >>>>
 >>>> And then imageB can be warped onto the extent of this image simply 
with:
 >>>> $ gdalwarp imageB.tiff domainA.tiff
 >>>>
 >>>> domainA.tiff will then contain the data from imageB, exactly 
co-located
 >>>> with imageA.
 >>>>
 >>>> The problem is that this only works when imageA contains a 
geotransform,
 >>>> and not when it contains only GCPs.
 >>>>
 >>>>
 >>>> So my question is:
 >>>> Is there any workaround to warp data from one image onto another image
 >>>> which is geolocated by GCPs only?
 >>>>
 >>>>
 >>>>
 >>>> With the -to option of gdalwarp it is possible to pass some parameters
 >>>> to GDALCreateGenImgProjTransformer2(). In the documentation of this
 >>>> function I find the following option:
 >>>>
 >>>> "METHOD: may have a value which is one of GEOTRANSFORM, 
GCP_POLYNOMIAL,
 >>>> GCP_TPS, GEOLOC_ARRAY, RPC to force only one geolocation method to be
 >>>> considered on the source dataset."
 >>>> http://www.gdal.org/gdal__alg_8h.html#94cd172f78dbc41d6f407d662914f2e3
 >>>>
 >>>> But there is no corresponding option to force use of 
GCP_POLYNOMIAL for
 >>>> the *destination* dataset. I guess this is a bad sign?
 >>>>
 >>>>
 >>>> Thank you for any help!
 >>>> Knut-Frode
 >>>>
 >>>> _______________________________________________
 >>>> gdal-dev mailing list
 >>>> gdal-dev at lists.osgeo.org
 >>>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
 >>> _______________________________________________
 >>> gdal-dev mailing list
 >>> gdal-dev at lists.osgeo.org
 >>> http://lists.osgeo.org/mailman/listinfo/gdal-dev
 >> _______________________________________________
 >> gdal-dev mailing list
 >> gdal-dev at lists.osgeo.org
 >> http://lists.osgeo.org/mailman/listinfo/gdal-dev
 >
 >