[gdal-dev] Warping onto an image with only GCPs with Python
Brian Case
rush at winkey.org
Fri Jan 13 14:40:23 EST 2012
Anton,
you can get the extents of the modis swath file from the metadata.
you have to use the the geolocation arrays to create a grid from the
swath file. the further left and right from the center of the image the
more area the single pixel in the swath file covers. also modis is a
pushbroom sensor. every 10 lines in the 1k file (20 in the 500m and 40
in the 250m) overlap the previous and next set. this is called the
bowtie effect. as it is right now gdal cannot properly deal with this.
short term I would use swath2grid to create tiffs of your data.
see http://trac.osgeo.org/gdal/ticket/4193
http://trac.osgeo.org/gdal/ticket/4192
http://trac.osgeo.org/gdal/ticket/4188
if you have any additional information or issues please add them
Brian
On Fri, 2012-01-13 at 10:33 +0100, Anton Korosov wrote:
> 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
> >
> >
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
More information about the gdal-dev
mailing list