warp introduces new values (using NearestNeighbour and Mode)

[kempenep]

warp introduces new values (using NearestNeighbour and Mode)

I expected nearest neighbour and mode not to produce new pixel values. I tried with the gdalwarp command line tool, Python bindings and C++ API. All produced similar results (new values introduced). The raster file on which I tried is the SCL file from the S2 image (classification file produced by sen2cor):

S2A_MSIL2A_20190329T095031_N0211_R079_T33TWN_20190329T124344.SAFE

Steps to reproduce the problem.

Using GDAL command line tool

gdalwarp -t_srs epsg:32632 -r near -tr 200 200 T33TWN_20190329T095031_SCL_20m.jp2 gdalwarped_near.tif

Using Python bindings:

`` from osgeo import gdal import numpy as np

fn='T33TWN_20190329T095031_SCL_20m.jp2'

warpfn1='/tmp/pythonwarp_nn.tif' warpfn2='/tmp/pythonwarp_mode.tif' warpfn3='/tmp/pythonwarp_cubic.tif' warpfn4='/tmp/pythonwarp_bilinear.tif'

ds = gdal.Open(fn) band = ds.GetRasterBand(1)

define different options

opt1 = gdal.WarpOptions(xRes=200, yRes=200, dstSRS='epsg:3035', resampleAlg=gdal.GRA_NearestNeighbour) opt2 = gdal.WarpOptions(xRes=200, yRes=200, dstSRS='epsg:3035', resampleAlg=gdal.GRA_Mode)

Read the data into numpy array

npa=np.array(ds.GetRasterBand(1).ReadAsArray())

perform warp with different options

gdal.Warp(warpfn1, ds, options=opt1) gdal.Warp(warpfn2, ds, options=opt2)

open warped datasets

ds1 = gdal.Open(warpfn1) ds2 = gdal.Open(warpfn2)

Read the data into numpy array

np1=np.array(ds1.GetRasterBand(1).ReadAsArray()) np2=np.array(ds2.GetRasterBand(1).ReadAsArray())

Get unique pixel values

print(np.unique(npa)) print(np.unique(np1)) print(np.unique(np2)) `` scl.zip

output (first: unique values from original SCL raster dataset, second and third: unique values of warped raster dataset using NN and mode): [ 0 2 3 4 5 6 7 8 9 10 11] [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14] [ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14]

Operating system: Debian 9 (64 bit)

GDAL version and provenance: 2.1.2 from Debian 9 stretch/main amd64 Packages

Thanks,

Pieter.

[rouault]

Ah this is a subtle issue... This is related to how JPEG2000 compression works, that you are subsampling (changing from 20m to 200m) and your raster is a classified one.

So this JPEG2000 image has 4 decomposition levels overviews), and as you are subsampling, GDAL will try to use the overviews whose spatial resolution is the closest to the target one you need (to save computation). However in your case, it turns out that this is a classification, and not imagery. But due how to wavelet compression works, JPEG2000 decomposition levels can have different pixel values that are not present in the full resolution image. So you have to instruct gdal.Warp() not to use those decomposition level, but just the full resolution image

So add overviewLevel = 'NONE' as an option of gdal.WarpOptions()

[kempenep]

Thanks for this clear answer. I guess the overviewLevel option is not yet available in GDAL 2.1.2. From gdalwarp_lib.cpp I understand the corresponding option in the C++ API would be nOvLevel (also not available in GDAL 2.1.2).

[rouault]

Ah, indeed overviewLevel has been added in GDAL 2.5 You can still use the command line style options gdal.Warp(out, in, options="-tr 200 200 -t_srs EPSG:3035 -r nearest -ovr NONE")

[kempenep]

Perfect, that works (for those interested re-using this line, use -r near instead of -r nearest). Thanks.

Op di 30 apr. 2019 om 12:25 schreef Even Rouault notifications@github.com:

Ah, indeed overviewLevel has been added in GDAL 2.5 You can still use the command line style options gdal.Warp(out, in, options="-tr 200 200 -t_srs EPSG:3035 -r nearest -ovr NONE")

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/OSGeo/gdal/issues/1484#issuecomment-487901439, or mute the thread https://github.com/notifications/unsubscribe-auth/AC3INSL6J2QU5WR5HEQI3U3PTANBHANCNFSM4HJBCEYA .