[PATCH] gdalwarp doesn't use overviews

[crschmidt]

It appears that gdalwarp does not take advantage of overviews, always choosing to read the entire source data when doing warping.

For the WMS driver, this means that warping any section of the map results in it reading the entire dataset at the lowest level from the dataset. In the case of OAM, this actually results in entirely black blocks being read for most areas, because the data isn't high enough resolution to actually have any pixels displayed per tile/block.

The end result is that gdalwarp against a WMS service essentially doesn't work. (Instead, one must project the coordinates, use gdal_translate, then gdalwarp the resulting image.)

To reproduce, do something like:

gdalwarp -te

233978.513603 899703.853477 237679.486397 902022.146523 -ts 803 503 -t_srs EPSG:2805 gdal_wms.xml foo.jpg

Using  http://openaerialmap.org/static/gdal_wms.xml

[bishop]

The overviews created by VRTWarpedDataset::IBuildOverviews also read entire the file. It's very slow in any case: if warped image have overvies or if they are created by IBuildOverviews

[bishop]

I tested this patch on tiff sample from  ftp://ftp.remotesensing.org/pub/geotiff/samples/misc/

[bishop]

This patch adds overview support for warped rasters with overviews (updated)

[bishop]

I have rewrote my patch for current trunk.

I added a couple of config options: GDAL_USE_SOURCE_OVERVIEWS and GDAL_SOURCE_OVERVIEWS_APROX.

If GDAL_USE_SOURCE_OVERVIEWS is set to ON (default is OFF) the Warped VRT Dataset will try to use source Raster overviews. If GDAL_USE_SOURCE_OVERVIEWS is set to OFF, the default GDAL behaviour is present.

Also I changed build scripts and dox file.

Need the advice, can I commit changes to trunk?

[rouault]

I didn't test the patch, but from review, I have a few remarks :

• GDAL_SOURCE_OVERVIEWS_APROX should be named GDAL_SOURCE_OVERVIEWS_APPROX (2 p's)
• I see indentation oddities. Make sure you use 4 spaces instead of tabs
• In VRTWarpedOverviewDataset::VRTWarpedOverviewDataset(), the calls to SetProjection?() and SetGCPs() don't need the VRTWarpedOverviewDataset:: qualifier
• In VRTWarpedDataset::IBuildOverviews(), at the end of the // Cleanup overviews section, before calling eErr = LoadSourceOverviews?(), I think that you need to add "papoOverviews = NULL; nOverviewCount = 0;" otherwise LoadSourceOverviews?() will work on a free'd pointer.
• I'm surprised that you can simply do SetBand?(i + 1, pOverviewDS->GetRasterBand?(i + 1)); in VRTWarpedOverviewDataset::VRTWarpedOverviewDataset(). The band has now 2 parents, the pOverviewDS and VRTWarpedOverviewDataset. So I would expect that when the VRTWarpedOverviewDataset is destroyed, the base destructor in GDALDataset will destroy the bands, which also belong to pOverviewDS and thus will be destroyed when pOverviewDS is destroyed... So my guess is that VRTWarpedOverviewDataset is actually never destroyed. (or that you're lucky to not observe crashes)
• In VRTWarpedOverviewDataset::VRTWarpedOverviewDataset(), I'm surprised that you can affect the georeferencing information (particularly the geotransform and GCP) of the pParentDS to the VRTWarpedOverviewDataset object, which has the same raster dimensions as the pOverviewDS, but not the dimensions of pParentDS. I may have not understood something, but I find that suspicious...
• Related to previous questions, why not directly providing pOverview to GDALAutoCreateWarpedVRT() ? What is the reason for VRTWarpedOverviewDataset ?
• Is the new code triggered by the existing Python autotests ? If not, it might be appropriate to extend them.

[bishop]

Replying to rouault:

I didn't test the patch, but from review, I have a few remarks : * GDAL_SOURCE_OVERVIEWS_APROX should be named GDAL_SOURCE_OVERVIEWS_APPROX (2 p's)

Fixed

* I see indentation oddities. Make sure you use 4 spaces instead of tabs

Fixed

* In VRTWarpedOverviewDataset::VRTWarpedOverviewDataset(), the calls to SetProjection?() and SetGCPs() don't need the VRTWarpedOverviewDataset:: qualifier

Fixed

* In VRTWarpedDataset::IBuildOverviews(), at the end of the // Cleanup overviews section, before calling eErr = LoadSourceOverviews?(), I think that you need to add "papoOverviews = NULL; nOverviewCount = 0;" otherwise LoadSourceOverviews?() will work on a free'd pointer.

Fixed

* I'm surprised that you can simply do SetBand?(i + 1, pOverviewDS->GetRasterBand?(i + 1)); in VRTWarpedOverviewDataset::VRTWarpedOverviewDataset(). The band has now 2 parents, the pOverviewDS and VRTWarpedOverviewDataset. So I would expect that when the VRTWarpedOverviewDataset is destroyed, the base destructor in GDALDataset will destroy the bands, which also belong to pOverviewDS and thus will be destroyed when pOverviewDS is destroyed... So my guess is that VRTWarpedOverviewDataset is actually never destroyed. (or that you're lucky to not observe crashes)

Fixed

* In VRTWarpedOverviewDataset::VRTWarpedOverviewDataset(), I'm surprised that you can affect the georeferencing information (particularly the geotransform and GCP) of the pParentDS to the VRTWarpedOverviewDataset object, which has the same raster dimensions as the pOverviewDS, but not the dimensions of pParentDS. I may have not understood something, but I find that suspicious...

As I know, this data (geotransform and GCP) use only for warping overview. While RasterIO GDAL use warp option of Dataset, not overview and everything works well. Maybe I mistaken, but in QGIS and wxGIS this works excelent on test images.

* Related to previous questions, why not directly providing pOverview to GDALAutoCreateWarpedVRT() ? What is the reason for VRTWarpedOverviewDataset ?

Unfortunatly pOverview has no Projection, geotransform and GCP's. I need some proxy dataset, which provides this information in GDALAutoCreateWarpedVRT while creation VRTWarpedDataset for overview. I try to use GDALProxyPoolDataset, but it more complicated and I cannot find how to SetBand?.

* Is the new code triggered by the existing Python autotests ? If not, it might be appropriate to extend them.

Please point me how and where I can check if Python autotests is triggered by code and where I should to add new tests for this code.

I rewrite the patch with fixed one.

[rouault]

The indentation changes in vrtdataset.cpp, vrtderivedrasterband.cpp, vrtdriver.cpp and vrtsources should be reverted. They add unrelated "noise" to the commit (and they should likely be fixed manually since the automatic fixing looks wrong). In order to keep the commit history clean, I find it better to separate indentation fix commits from semantic change commits for old code. For new code, it is better to do it right at first time.

I should likely investigate in-depth the issue with the georeferencing info on VRTWarpedOverviewDataset. I still think it is wrong. For example, the pixel size components of a geotransform depend of the raster dimensions. You cannot just copy a geotransform of a full resolution dataset to its overviews (otherwise they will be misplaced if someone tries to access the geotransform on the overviews). So if the end result is correct, this might be the evidence that some code requires things that it doesn't really use. I also think that the better fix would be to make sure that the overview dataset has valid geotransform, GCP, projection, which would make VRTWarpedOverviewDataset unnecesserary, but that might be more involved to track (the changes would be in the GTIFF, HFA and potentially other drivers that generate a dataset overviews).

As far as checking if the code is triggered, one way is to follow the "Test coverage" section of http://trac.osgeo.org/gdal/wiki/AutotestStatus. This is what I use to generate the coverage reports at  http://even.rouault.free.fr.

A simpler method is to add some fprintf(log, "....") statements in your code and see if they are run. As far as extending the autotest suite, this could likely go in autotest/gdrivers/vrtwarp.py (but they are other scripts that likely trigger warped VRTs, such as autotest/warp/warp.py).

[bishop]

Replying to rouault:

The indentation changes in vrtdataset.cpp, vrtderivedrasterband.cpp, vrtdriver.cpp and vrtsources should be reverted. They add unrelated "noise" to the commit (and they should likely be fixed manually since the automatic fixing looks wrong). In order to keep the commit history clean, I find it better to separate indentation fix commits from semantic change commits for old code. For new code, it is better to do it right at first time.

Fixed

I should likely investigate in-depth the issue with the georeferencing info on VRTWarpedOverviewDataset. I still think it is wrong. For example, the pixel size components of a geotransform depend of the raster dimensions. You cannot just copy a geotransform of a full resolution dataset to its overviews (otherwise they will be misplaced if someone tries to access the geotransform on the overviews). So if the end result is correct, this might be the evidence that some code requires things that it doesn't really use. I also think that the better fix would be to make sure that the overview dataset has valid geotransform, GCP, projection, which would make VRTWarpedOverviewDataset unnecesserary, but that might be more involved to track (the changes would be in the GTIFF, HFA and potentially other drivers that generate a dataset overviews).

I think that changing overviews behaviour in most formats is not good idea. I can explain how to recompute GeoTransformation?, but GCP's cannot be recomputed (i.e. the GCP's on neighbour pixels of main raster, in overview will point on same pixel). Also we will have some extra information which we need keep in mind while the main raster modifying. I think for such unique situation (warped dataset) the kind of hack is suitable. Also this is not default behavior of warped dataset according to config key GDAL_USE_SOURCE_OVERVIEWS. So, which developers who decide to use new functionality will do it on their own risk.

As far as checking if the code is triggered, one way is to follow the "Test coverage" section of http://trac.osgeo.org/gdal/wiki/AutotestStatus. This is what I use to generate the coverage reports at  http://even.rouault.free.fr. A simpler method is to add some fprintf(log, "....") statements in your code and see if they are run. As far as extending the autotest suite, this could likely go in autotest/gdrivers/vrtwarp.py (but they are other scripts that likely trigger warped VRTs, such as autotest/warp/warp.py).

Would you mind to run the tests for me? I don't have an Linux machine at hand. I'm afraid coverage.patch will not work on Windows...

[rouault]

Hum, after some preliminary testing, there are various issues :

1) Minor, but in the provided patch the GDAL_USE_SOURCE_OVERVIEWS is turned ON, whereas it is documented to be OFF. Actually, once the patch has settled down and been made rock solid, I think turning it ON might be appropriate.

2) The pSrcDS->BuildOverviews?() in VRTWarpedDataset::IBuildOverviews() is not appropriate. People who use VRT use it because it is fast. But currently BuildOverviews?() is called at dataset opening if a <OverviewList?> is defined in the VRT XML. And calling BuildOverviews?() on a large geotiff is very slow if the overview already exist. I've noticed that on a 10474x4951 geotiff with 7 overview levels.

3) Disabling that call speeds up things a lot, but unfortunately generates many Valgrind errors and finish by segfaulting. The first error reported is :

==31429== Invalid read of size 4
==31429==    at 0x57FFB7A: GDALDataset::GetRasterXSize() (gdaldataset.cpp:583)
==31429==    by 0x57DD2D3: VRTWarpedOverviewDataset::VRTWarpedOverviewDataset(GDALDataset*, GDALDataset*) (vrtwarpedoverview.cpp:42)
==31429==    by 0x57DACDC: VRTWarpedDataset::LoadSourceOverviews() (vrtwarped.cpp:446)
==31429==    by 0x57DB583: VRTWarpedDataset::IBuildOverviews(char const*, int, int*, int, int*, int (*)(double, char const*, void*), void*) (vrtwarped.cpp:744)
==31429==    by 0x5800776: GDALDataset::BuildOverviews(char const*, int, int*, int, int*, int (*)(double, char const*, void*), void*) (gdaldataset.cpp:1396)
==31429==    by 0x57DBF3D: VRTWarpedDataset::XMLInit(CPLXMLNode*, char const*) (vrtwarped.cpp:1029)
==31429==    by 0x57C6DC9: VRTDataset::OpenXML(char const*, char const*, GDALAccess) (vrtdataset.cpp:836)
==31429==    by 0x57C6AFF: VRTDataset::Open(GDALOpenInfo*) (vrtdataset.cpp:757)
==31429==    by 0x5801DD6: GDALOpenInternal(GDALOpenInfo&, char const* const*) (gdaldataset.cpp:2257)
==31429==    by 0x5801BBC: GDALOpenInternal(char const*, GDALAccess, char const* const*) (gdaldataset.cpp:2215)
==31429==    by 0x5801B76: GDALOpen (gdaldataset.cpp:2206)
==31429==    by 0x403C81: main (gdalinfo.c:178)
==31429==  Address 0x172dea94 is 52 bytes inside a block of size 264 free'd
==31429==    at 0x4C283A4: operator delete(void*) (vg_replace_malloc.c:480)
==31429==    by 0x57DD762: VRTWarpedOverviewDataset::~VRTWarpedOverviewDataset() (vrtwarpedoverview.cpp:74)
==31429==    by 0x580240B: GDALClose (gdaldataset.cpp:2448)
==31429==    by 0x57DA97A: VRTWarpedDataset::CloseDependentDatasets() (vrtwarped.cpp:352)
==31429==    by 0x57DA7E8: VRTWarpedDataset::~VRTWarpedDataset() (vrtwarped.cpp:298)
==31429==    by 0x580240B: GDALClose (gdaldataset.cpp:2448)
==31429==    by 0x57DB520: VRTWarpedDataset::IBuildOverviews(char const*, int, int*, int, int*, int (*)(double, char const*, void*), void*) (vrtwarped.cpp:735)
==31429==    by 0x5800776: GDALDataset::BuildOverviews(char const*, int, int*, int, int*, int (*)(double, char const*, void*), void*) (gdaldataset.cpp:1396)
==31429==    by 0x57DBF3D: VRTWarpedDataset::XMLInit(CPLXMLNode*, char const*) (vrtwarped.cpp:1029)
==31429==    by 0x57C6DC9: VRTDataset::OpenXML(char const*, char const*, GDALAccess) (vrtdataset.cpp:836)
==31429==    by 0x57C6AFF: VRTDataset::Open(GDALOpenInfo*) (vrtdataset.cpp:757)
==31429==    by 0x5801DD6: GDALOpenInternal(GDALOpenInfo&, char const* const*) (gdaldataset.cpp:2257)

[bishop]

Replying to rouault:

Hum, after some preliminary testing, there are various issues : 1) Minor, but in the provided patch the GDAL_USE_SOURCE_OVERVIEWS is turned ON, whereas it is documented to be OFF. Actually, once the patch has settled down and been made rock solid, I think turning it ON might be appropriate.

Pardon me, I set ON for testing. Fixed

2) The pSrcDS->BuildOverviews?() in VRTWarpedDataset::IBuildOverviews() is not appropriate. People who use VRT use it because it is fast. But currently BuildOverviews?() is called at dataset opening if a <OverviewList?> is defined in the VRT XML. And calling BuildOverviews?() on a large geotiff is very slow if the overview already exist. I've noticed that on a 10474x4951 geotiff with 7 overview levels.

The idea to create overviews for source Dataset came from GeoTransformation? with 45 degree rotation. In such cases the output overviews was rather big and taken a lot of time to generate. Building overviews for Source Dataset was faster. This is mainly need for programatically warped datasets. So, I decide to drop this controversial feature. Also, there is another way to create overviews for source dataset.

The main aim of this patch was directed to programatically use of warped datasets. It's very good, that you tested another using of them.

[bishop]

This patch adds overview support for warped rasters with overviews (rewrite for current trunk)

[rouault]

I'm seeing misplacement issues when displaying in QGIS.

How to replicate:

$ gdalinfo world_4326.tif
Driver: GTiff/GeoTIFF
Files: world_4326.tif
Size is 10474, 4951
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-179.999999974943762,85.084059047963109)
Pixel Size = (0.034370672254824,-0.034370672254824)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  85.0840590) (180d 0' 0.00"W, 85d 5' 2.61"N)
Lower Left  (-180.0000000, -85.0851393) (180d 0' 0.00"W, 85d 5' 6.50"S)
Upper Right ( 179.9984212,  85.0840590) (179d59'54.32"E, 85d 5' 2.61"N)
Lower Right ( 179.9984212, -85.0851393) (179d59'54.32"E, 85d 5' 6.50"S)
Center      (  -0.0007894,  -0.0005401) (  0d 0' 2.84"W,  0d 0' 1.94"S)
Band 1 Block=10474x1 Type=Byte, ColorInterp=Red
Band 2 Block=10474x1 Type=Byte, ColorInterp=Green
Band 3 Block=10474x1 Type=Byte, ColorInterp=Blue

$ gdaladdo world_4326.tif 2 4 8 16 -ro

$ GDAL_USE_SOURCE_OVERVIEWS=YES gdalwarp -of vrt world_4326.tif test.vrt -t_srs EPSG:3857 -overwrite

$ GDAL_USE_SOURCE_OVERVIEWS=YES qgis test.vrt

Then if you center the map on Paris at full resolution and progressively zoom out, the center at some point jumps to somewhere in Africa. The same procedure without GDAL_USE_SOURCE_OVERVIEWS=YES doesn't show this (although being much slower).

Other point: while analyzing call stacks in the debugger, I see that the complete removal of the special case for GDAL_USE_SOURCE_OVERVIEWS=YES is perhaps a bit extreme since now in VRTWarpedDataset::XMLInit(), the overview datasets that are built are the "slow" VRTWarpedDataset based on the full resolution source dataset, but those overviews have themselves overviews based on VRTWarpedOverviewDataset... So I assume that we don't take advantage of the full speed-up that could be taken advantage of.

[bishop]

Thank's Even for good example. I reproduced the bug. Unfortunatly I need much more time to fix it. The other errors will also be fixed, but they connected with the main problem.