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

[Brian Case]

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