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",50,
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.000,4390830.000)
Pixel Size = (0.150,-0.150)
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
retur
