[gdal-dev] Issue with resampling algorithm in gdal_array.DatasetReadAsArray python wrapper function
Even Rouault
even.rouault at spatialys.com
Thu Feb 10 08:16:23 PST 2022
Ok, so yes the issue is likely a bug in the ECW driver. Quickly looking
at the code I do see something potentially suspicious. Please file a
ticket to https://github.com/OSGeo/gdal/issues/new with all the below
information.
Le 10/02/2022 à 17:01, Miguel A. Manso a écrit :
> 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
>>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
More information about the gdal-dev
mailing list