[gdal-dev] Interpolation when using the WMS format

[Peter Hopfgartner]

On 09/18/2014 10:30 PM, Even Rouault wrote:
> Le jeudi 18 septembre 2014 22:23:06, Peter Hopfgartner a écrit :
>> Salut Evan,
>>
>> thanks for you recipe. Unfortunately, what I am really interested in is to
>> use this from MapServer. Could you give me a pointer in the source where
>> the resampling takes place?
> In gcore/rasterio.cpp : GDALRasterBand::IRasterIO(). I must warn you : quite
> tricky code
>
> But perhaps you don't need code change. In MapServer did you try to define
> PROCESSING "RESAMPLE=BILINEAR" ?
>
>> Thanks,
>>
>> Peter
>>
>> Am 18.09.2014 um 22:11 schrieb Even Rouault <even.rouault at spatialys.com>:
>>> Le jeudi 18 septembre 2014 17:40:11, Peter Hopfgartner a écrit :
>>>> Hi, everybody
>>>>
>>>> we are using the GDAL WMS driver for getting maps from a public WMTS
>>>> service, http://www.basemap.at. Anyway, the resulting image does not
>>>> look smooth. In particular, I would guess that the the image was scaled
>>>> with something like a nearest neighbor algorithm.
>>>> So my question is simply, if I can switch the interpolation method in
>>>> some place.
>>>>
>>>> The XML file, taken from
>>>> http://gisforge.wordpress.com/2012/11/06/basemap-at-wmts-layer-in-qgis/,
>>>> is:
>>>>
>>>> <GDAL_WMS>
>>>>
>>>>   <Service name="TMS">
>>>>
>>>> <ServerUrl>http://maps.wien.gv.at/basemap/geolandbasemap/normal/google38
>>>> 57/ ${z}/${y}/${x}.jpeg</ServerUrl> </Service>
>>>>
>>>>   <DataWindow>
>>>>   
>>>>    <UpperLeftX>-20037508.34</UpperLeftX>
>>>>    <UpperLeftY>20037508.34</UpperLeftY>
>>>>    <LowerRightX>20037508.34</LowerRightX>
>>>>    <LowerRightY>-20037508.34</LowerRightY>
>>>>    <TileLevel>18</TileLevel>
>>>>    <TileCountX>1</TileCountX>
>>>>    <TileCountY>1</TileCountY>
>>>>    <YOrigin>top</YOrigin>
>>>>   
>>>>   </DataWindow>
>>>>   <Projection>EPSG:3857</Projection>
>>>>   <BlockSizeX>256</BlockSizeX>
>>>>   <BlockSizeY>256</BlockSizeY>
>>>>   <BandsCount>3</BandsCount>
>>>>   <Cache />
>>>>
>>>> </GDAL_WMS>
>>>>
>>>> The command is:
>>>>
>>>> gdal_translate -srcwin 35500000 23500000 100000 100000 -of JPEG -outsize
>>>> 250 250 basemap_at.xml basemap_at.jpg
>>> Peter,
>>>
>>> Resampling done by gdal_translate is nearest neighbour, so not always
>>> that great.
>>>
>>> What you can do is select an outsize that is compatible with the source
>>> overviews. In TMS, the ratio between sucessive overviews work with factor
>>> of 2, so you should select an outsize of 10000 / 2^N. The closest is
>>> 10000 / 256 = 390.625 -> 391
>>>
>>> So
>>>
>>> gdal_translate -srcwin 35500000 23500000 100000 100000 -outsize
>>> 391 391 basemap_at.xml tmp.tif
>>>
>>> and then:
>>>
>>> gdalwarp tmp.tif tmp2.tif -r lanczos -ts 250 250 (you can try different
>>> resampling kernel)
>>> gdal_translate tmp2.tif basemap_at.jpg -of JPEG
>>>
>>> Best regards,
>>>
>>> Even
I've tried to simulate the result of different algorithms to further 
investigate the idea of having a configurable resampling step. The 
original tile is 256 x 256 pixels and I've rescaled it to rescale it to 
200x200 using Gimp.

The original tile was taken from 
http://maps.wien.gv.at/basemap/geolandbasemap/normal/google3857/8/89/137.jpg, 
a copy is stored at 
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.jpg.

Results
no interpolation, (nearest neighbor?): 
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.resampled_to_200.px-noint.png
bilinear: 
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.resampled_to_200.px-linear.png
bicubic: 
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.resampled_to_200.px-cubic.png
lanczos: 
https://dl.dropboxusercontent.com/u/4935329/resampling_tests/137.resampled_to_200.px-lanczos.png

I'd say that the quality of the resampling respects the obove order, 
with lanczos being the best.

Regards,

Peter

-- 
Peter Hopfgartner
R3 GIS
web  : http://www.r3-gis.com