After looking into GDAL code the overview logic seems similar to that in Mapnik. I was thinking there needs to be some recursion to go through nested datasets, but I didn't find such.

It would really help if you can make a simple test case which I can run locally.

I agree that the way how Mapnik GDAL plugin handles overviews is not how some users would expect.

I've created an example. You can find it at https://gist.github.com/cpaulik/34853bebc07174edc78233be51fc8bb6

Thanks for the test case, I can reproduce it.

Mapnik can see only one overview, as well as gdalinfo. Other overviews are hidden behind abstraction of another data set.

$ gdalinfo global.vrt
Driver: VRT/Virtual Raster
Files: global.vrt
       global.vrt.ovr
       5x5.vrt
Size is 6000, 6000
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 = (-10.000000000000000,40.000000000000000)
Pixel Size = (0.001000000000000,-0.001000000000000)
Corner Coordinates:
Upper Left  ( -10.0000000,  40.0000000) ( 10d 0' 0.00"W, 40d 0' 0.00"N)
Lower Left  ( -10.0000000,  34.0000000) ( 10d 0' 0.00"W, 34d 0' 0.00"N)
Upper Right (  -4.0000000,  40.0000000) (  4d 0' 0.00"W, 40d 0' 0.00"N)
Lower Right (  -4.0000000,  34.0000000) (  4d 0' 0.00"W, 34d 0' 0.00"N)
Center      (  -7.0000000,  37.0000000) (  7d 0' 0.00"W, 37d 0' 0.00"N)
Band 1 Block=128x128 Type=Float64, ColorInterp=Gray
  Overviews: 750x750

Mapnik loops through all accessible overviews (just one in this case). Finds no suitable overview and instead of requesting raster in resolution of the map (256 x 256), requests full resolution of the raster for given extent.

The following python snippet is equivalent to the RasterIO called from Mapnik:

from osgeo import gdal
ds = gdal.Open('global.vrt')
band = ds.GetRasterBand(1)
data = band.ReadRaster(xoff=0, yoff=0,
                       xsize=257*4, ysize=257*4,
                       buf_xsize=257*4, buf_ysize=257*4,
                       buf_type=gdal.GDT_Float64)

Just a note: Mapnik adds 4 pixels margin, which is the reason why there is 257 * 4 instead of 256 * 4.

This is really bad for your use case. What options do we have? Make this behaviour optional?

What do you think, @flippmoke ?

I would have to look at the code for why there is a margin like this but my assumption would be that it is a slightly bigger buffer so that we can do better resampling if necessary on other code paths. I am not certain here at all though, it could simply be a bug. Do we know what line of code this comes from?

@flippmoke No, no, the margin is not the problem here at all. Margin is correct as it is. I just tried to explain those values.

The problem here is that Mapnik does not use overviews, instead requesting raster 16 times bigger than needed. There is no way how Mapnik can do that with current implementation as those overviews are hidden in nested VRT. But GDAL itself can use those overviews, the old implementation of Mapnik GDAL plugin as well.

Okay sorry, I was slightly confused here...

I put together a PR that might address this and gives me at least something visual to represent the problem. I have not tested it at all against this testcase and we likely should add a test case here.

Finally, we might want to consider if it is possible to make an endless loop here? Can you have a VRT overview that points to a VRT? This could be a security concern if that is the case.

@flippmoke I doubt that your PR addresses this problem since I think it is not possible to get at the nested overviews using a GetOverview call on the first dataset. It is that by using the correct size and buffer combination in RasterIO the nested gdal drivers will use overviews of all levels as appropriate.

I've heard that mapnik wants to do it's own resampling in some cases and this is the reason for implementing all this logic in mapnik. But as I said above. I think most users would expect the gdal plugin to work like gdal would.

Thanks for the pull request, @flippmoke, but it doesn't work, unfortunately. As @cpaulik suggests, we should not to try to outsmart GDAL and rather call RasterIO with expected parameters.

The current logic doesn't solve memory issues nor it solves resampling as overviews are already resampled out of Mapnik's control.

I am fine with doing resampling via GDAL, we will just need to have some sort of conversion between the different expected resampling types so that the correct one is chosen. Additionally, I am not sure how in the later steps if there are any concerns.

@rouault do you have any suggestions here?

@rouault do you have any suggestions here?

This is a bit vague as a question ;-) You can specify a resampling method to the extra argument of GDALRasterIOEx()

@flippmoke The resampling type is chosen when the overview is created so you have not choice there on mapniks side if I understand correctly.

If there are no overviews then it is a different story of course. But how would a mapnik user choose a resampling method? I don't see any parameter listed in the wiki https://github.com/mapnik/mapnik/wiki/GDAL

EDIT: Ah sorry, it seems to be in the RasterSymbolyzer https://github.com/mapnik/mapnik/wiki/RasterSymbolizer

@rouault is there a way to access an overview in a TIF from an overview in a VRT? Should we just trust GDAL in this case? Would it be more wise for us to simply have GDAL do all the band selection and resampling?

is there a way to access an overview in a TIF from an overview in a VRT?

Not directly since VRT abstracts away those details, but if you issue a downsampling RasterIO() request on the VRT, it will forward it to the TIF, and the TIF will then use its overviews when appropraite.

@flippmoke I've only found the three resampling types in the RasterSymbolizer documentation. Are those the ones we have to map to gdal resampling methods?

And we should probably add some kind of warning in the docs that this setting will be ignored if overviews exist.

@cpaulik all of the possible ways we resample https://github.com/mapnik/mapnik/blob/master/include/mapnik/image_scaling.hpp#L42-L61

That makes it more interesting 😄

That would be a first mapping without looking up all the resampling methods in detail.

Would you be ok in saying the resampling methods from mapnik that are not supported by gdal are not supported by the gdal plugin? Or would you prefer a different data access logic based on chosen resampling method?

What I mean is that if e.g. SCALING_BESSEL is chosen in mapnik then mapnik will request the full resolution dataset and then resample itself but if e.g. SCALING_NEAR is chosen then mapnik would call RasterIO with the final buffer size so that gdal resamples / or uses overviews if present.

@cpaulik just some history

The chart above is part of the reason we have not done resampling in GDAL in the past the expected result if you are using Mapnik resampling can be very different then GDAL. Therefore the choice was made to try to resample using Mapnik where possible and then search through the overviews to find the best overview to select the source data from for resampling. Part of the goal has been for the raster plugin and the GDAL plugin to reproduce the same results if possible.

I am worried that in solving one problem we could be creating another and do something that is another big shift in the GDAL plugin/mapnik. Should we do this?

@flippmoke Thank you for the background.

My thoughts on this are that the current solution does not solve the problem. If any resampling method is selected that does not produce the same result as gdal resampling then mapnik will produce inconsistent results based on zoom level. Especially if overviews are present. In addition it is very hard to find out why mapnik does not work as expected with a gdal datasource. I've spend several hours figuring out why mapnik 3 does not work properly with our VRT files.

There are several approaches to this problem that I can think of that would make it clearer from a user perspective:

1. The gdal plugin does exactly what gdal would do. So the users can test with gdal, or any other library/application that uses gdal and then switch to mapnik without a lot of surprises.
2. The gdal plugin is smart enough to use the gdal RasterIO differently based on the selected settings. So if a resampling method is selected that is not supported by gdal then the resampling is done in mapnik and overviews are ignored.
3. A combination of 1 and 2 where the user can select a e.g. pure-gdal mode or a mapnik mode.

Part of the goal has been for the raster plugin and the GDAL plugin to reproduce the same results if possible.

I think it'd be better for the GDAL plugin to delegate as much as possible to libgdal. It will be less surprising for users if it works similar to other gdal tools, rather than an unrelated raster plugin.

If a resampling filter unsupported by GDAL is selected, the plugin should issue a warning and select a supported filter with similar characteristics (e.g. when a Sinc-family filter is selected, fallback to Lanczos).

If there is actual demand for resampling a GDAL datasource on mapnik side, it should be opt-in via parameters. And then perhaps another parameter could provide information about nested overviews.

If a resampling filter unsupported by GDAL is selected, the plugin should issue a warning and select a supported filter with similar characteristics (e.g. when a Sinc-family filter is selected, fallback to Lanczos).

That's a good compromise solution and would work beautifully for our use case.

If I can be any help in implementing this please tell me. I don't know the mapnik code base very well so I don't feel comfortable with making a PR for this at the moment.

I don't know the mapnik code base very well so I don't feel comfortable with making a PR for this at the moment.

@cpaulik updating PRs until it works is how most everyone in this ticket learned the mapnik code base 😄

@flippmoke Ok, If you are patient with me then I can try.

I was trying to get started and I'm curious to know if there is some documentation about the mapnik internals and how the different classes work with each other. There also does not seem to be any documentation in the code about what all the parameters mean. Am I missing something or does it really not exist?

Am I missing something or does it really not exist?

Unfortunately there really is no documentation about Mapnik's internals.