[gdal-dev] Issue with resampling algorithm in gdal_array.DatasetReadAsArray python wrapper function

[Miguel A. Manso]

Even

Thank you for your comments
I have created a new environment, installed the latest version (3.4.1) 
had version 3.0.4 and the problem persisted.
I have done another test which may give the clue as to where the problem 
lies. My datasource is an ECW file and if I use this datasource to 
extract an area with the gdal_array.DatasetReadAsArray function telling 
it to use resample methods other than NearestNeighbour it does not do it 
right.
I have cut a slice and stored it as a geotiff, and using this datasource 
for the same test has already worked the resample algorithms. It works 
with both versions of the libraries.
I'm sorry I can't extract a piece of image to share it as I have no way 
to store it in ecw format.

The code used:


ds = gdal.Open("PNOA_MA_OF_ETRS89_HU31_h50_0699.ecw")

x = gdal_array.DatasetReadAsArray( ds, 1, 1, 256, 256, buf_obj=None, 
buf_xsize= 384, buf_ysize=384, buf_type=None, 
resample_alg=gdal.GRIORA_Lanczos)

src_srs = ogr.osr.SpatialReference()
src_srs.ImportFromEPSG(25830)
wkt_projection = src_srs.ExportToWkt()

driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create( "recorte.tif", 384, 384, 3, gdal.GDT_Byte )
dataset.SetGeoTransform(( x_min, 0.1, 0, y_max, 0, -0.1))
dataset.SetProjection(wkt_projection)
bands = array.shape[0]
for band in range(bands):
         dataset.GetRasterBand(band+1).WriteArray( array[ band,:, : ] )
dataset.FlushCache()
bands = None
dataset = None
driver = None


gdalinfo over ECW:

Driver: ECW/ERDAS Compressed Wavelets (SDK 5.3)
Files: PNOA_MA_OF_ETRS89_HU31_h50_0699.ecw
Size is 191800, 124133
Coordinate System is:
PROJCRS["NUTM31",
     BASEGEOGCRS["ETRS89",
         DATUM["European Terrestrial Reference System 1989",
             ELLIPSOID["GRS 1980",6378137,298.257222101,
                 LENGTHUNIT["metre",1]]],
         PRIMEM["Greenwich",0,
             ANGLEUNIT["degree",0.0174532925199433]],
         ID["EPSG",4258]],
     CONVERSION["UTM zone 31N",
         METHOD["Transverse Mercator",
             ID["EPSG",9807]],
         PARAMETER["Latitude of natural origin",0,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8801]],
         PARAMETER["Longitude of natural origin",3,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8802]],
         PARAMETER["Scale factor at natural origin",0.9996,
             SCALEUNIT["unity",1],
             ID["EPSG",8805]],
         PARAMETER["False easting",500000,
             LENGTHUNIT["Metre",1],
             ID["EPSG",8806]],
         PARAMETER["False northing",0,
             LENGTHUNIT["Metre",1],
             ID["EPSG",8807]],
         ID["EPSG",16031]],
     CS[Cartesian,2],
         AXIS["easting",east,
             ORDER[1],
             LENGTHUNIT["Metre",1]],
         AXIS["northing",north,
             ORDER[2],
             LENGTHUNIT["Metre",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (483820.000000000000000,4390830.000000000000000)
Pixel Size = (0.150000000000000,-0.150000000000000)
Metadata:
   COLORSPACE=RGB
   COMPRESSION_RATE_TARGET=6
   VERSION=2
Corner Coordinates:
Upper Left  (  483820.000, 4390830.000) (  2d48'40.90"E, 39d40' 1.67"N)
Lower Left  (  483820.000, 4372210.050) (  2d48'42.54"E, 39d29'57.69"N)
Upper Right (  512590.000, 4390830.000) (  3d 8'48.42"E, 39d40' 1.89"N)
Lower Right (  512590.000, 4372210.050) (  3d 8'47.15"E, 39d29'57.90"N)
Center      (  498205.000, 4381520.025) (  2d58'44.75"E, 39d35' 0.22"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
   Description = Red
   Overviews: 95900x62066, 47950x31033, 23975x15516, 11987x7758, 
5993x3879, 2996x1939, 1498x969, 749x484, 374x242
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
   Description = Green
   Overviews: 95900x62066, 47950x31033, 23975x15516, 11987x7758, 
5993x3879, 2996x1939, 1498x969, 749x484, 374x242
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
   Description = Blue
   Overviews: 95900x62066, 47950x31033, 23975x15516, 11987x7758, 
5993x3879, 2996x1939, 1498x969, 749x484, 374x242


Thank you

Miguel A


El 10/02/2022 a las 14:06, Even Rouault escribió:
> Miguel,
>>
>> I am trying to use the gdal_array.DatasetReadAsArray function to 
>> extract a window from a 3B gtif file by indicating the dataset, the 
>> upper left corner image coordinates, the window size (256x256) at 
>> actual scale in the dataset, the desired buffer size (larger than 
>> actual size 341x341) and indicating the resample algorithm.
>> It only works if I select the NearestNeighbour algorithm. In the rest 
>> of the cases (Bilinear, bicubic, Lanczos, etc...) what the function 
>> returns me is a mosaic that repeats 9 times (3x3) the original image 
>> in reduced size.
>>
>> What am I doing wrong or what have I interpreted incorrectly in the API?
>
> First you should try with the latest release (3.4.1). If you can still 
> reproduce the issue, then strip down your code to a minimum standalone 
> reproducer script (using only GDAL & numpy) and with a minimum 
> reproducing dataset, and provide it here.
>
> Even
>
>>
>> Thanks & Best regards
>>
>> Miguel A. Manso
>>
>> MERCATOR Research Group
>>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>