922 строки
30 KiB
Python
Исполняемый файл
922 строки
30 KiB
Python
Исполняемый файл
#!/usr/bin/env pytest
|
|
# -*- coding: utf-8 -*-
|
|
###############################################################################
|
|
# $Id$
|
|
#
|
|
# Project: GDAL/OGR Test Suite
|
|
# Purpose: Test core numeric operations and statistics calculations
|
|
# Author: Mateusz Loskot <mateusz@loskot.net>
|
|
#
|
|
###############################################################################
|
|
# Copyright (c) 2007, Mateusz Loskot <mateusz@loskot.net>
|
|
# Copyright (c) 2009-2012, Even Rouault <even dot rouault at spatialys.com>
|
|
#
|
|
# This library is free software; you can redistribute it and/or
|
|
# modify it under the terms of the GNU Library General Public
|
|
# License as published by the Free Software Foundation; either
|
|
# version 2 of the License, or (at your option) any later version.
|
|
#
|
|
# This library is distributed in the hope that it will be useful,
|
|
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
# Library General Public License for more details.
|
|
#
|
|
# You should have received a copy of the GNU Library General Public
|
|
# License along with this library; if not, write to the
|
|
# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
|
|
# Boston, MA 02111-1307, USA.
|
|
###############################################################################
|
|
|
|
import math
|
|
import os
|
|
import shutil
|
|
import struct
|
|
|
|
import gdaltest
|
|
import pytest
|
|
|
|
from osgeo import gdal
|
|
|
|
###############################################################################
|
|
# Test handling NaN with GDT_Float32 data
|
|
|
|
|
|
def test_stats_nan_1():
|
|
|
|
gdaltest.gtiff_drv = gdal.GetDriverByName("GTiff")
|
|
if gdaltest.gtiff_drv is None:
|
|
pytest.skip()
|
|
|
|
stats = (50.0, 58.0, 54.0, 2.5819888974716)
|
|
|
|
shutil.copyfile("data/nan32.tif", "tmp/nan32.tif")
|
|
|
|
t = gdaltest.GDALTest("GTiff", "tmp/nan32.tif", 1, 874, filename_absolute=1)
|
|
ret = t.testOpen(check_approx_stat=stats, check_stat=stats)
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("tmp/nan32.tif")
|
|
|
|
return ret
|
|
|
|
|
|
###############################################################################
|
|
# Test handling NaN with GDT_Float64 data
|
|
|
|
|
|
def test_stats_nan_2():
|
|
|
|
if gdaltest.gtiff_drv is None:
|
|
pytest.skip()
|
|
|
|
stats = (50.0, 58.0, 54.0, 2.5819888974716)
|
|
|
|
shutil.copyfile("data/nan64.tif", "tmp/nan64.tif")
|
|
|
|
t = gdaltest.GDALTest("GTiff", "tmp/nan64.tif", 1, 4414, filename_absolute=1)
|
|
ret = t.testOpen(check_approx_stat=stats, check_stat=stats)
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("tmp/nan64.tif")
|
|
|
|
return ret
|
|
|
|
|
|
###############################################################################
|
|
# Test stats on signed byte (#3151)
|
|
|
|
|
|
def test_stats_signedbyte():
|
|
|
|
if gdaltest.gtiff_drv is None:
|
|
pytest.skip()
|
|
|
|
stats = (-128.0, 127.0, -0.2, 80.64)
|
|
|
|
shutil.copyfile("data/stats_signed_byte.img", "tmp/stats_signed_byte.img")
|
|
|
|
t = gdaltest.GDALTest(
|
|
"HFA", "tmp/stats_signed_byte.img", 1, 11, filename_absolute=1
|
|
)
|
|
ret = t.testOpen(check_approx_stat=stats, check_stat=stats, skip_checksum=1)
|
|
|
|
gdal.GetDriverByName("HFA").Delete("tmp/stats_signed_byte.img")
|
|
|
|
return ret
|
|
|
|
|
|
###############################################################################
|
|
# Test return of GetStatistics() when we don't have stats and don't
|
|
# force their computation (#3572)
|
|
|
|
|
|
def test_stats_dont_force():
|
|
|
|
if os.path.exists("data/byte.tif.aux.xml"):
|
|
gdal.Unlink("data/byte.tif.aux.xml")
|
|
ds = gdal.Open("data/byte.tif")
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 0)
|
|
assert stats == [0, 0, 0, -1], "did not get expected stats"
|
|
|
|
|
|
###############################################################################
|
|
# Test statistics when stored nodata value doesn't accurately match the nodata
|
|
# value used in the imagery (#3573)
|
|
|
|
|
|
def test_stats_approx_nodata():
|
|
|
|
shutil.copyfile("data/minfloat.tif", "tmp/minfloat.tif")
|
|
try:
|
|
os.remove("tmp/minfloat.tif.aux.xml")
|
|
except OSError:
|
|
pass
|
|
|
|
ds = gdal.Open("tmp/minfloat.tif")
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
md = ds.GetRasterBand(1).GetMetadata()
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
os.remove("tmp/minfloat.tif.aux.xml")
|
|
|
|
ds = gdal.Open("tmp/minfloat.tif")
|
|
minmax = ds.GetRasterBand(1).ComputeRasterMinMax()
|
|
ds = None
|
|
|
|
os.remove("tmp/minfloat.tif")
|
|
|
|
if nodata != -3.4028234663852886e38:
|
|
print("%.18g" % nodata)
|
|
pytest.fail("did not get expected nodata")
|
|
|
|
assert stats == [-3.0, 5.0, 1.0, 4.0], "did not get expected stats"
|
|
|
|
assert md == {
|
|
"STATISTICS_MEAN": "1",
|
|
"STATISTICS_MAXIMUM": "5",
|
|
"STATISTICS_MINIMUM": "-3",
|
|
"STATISTICS_STDDEV": "4",
|
|
"STATISTICS_VALID_PERCENT": "50",
|
|
}, "did not get expected metadata"
|
|
|
|
assert minmax == (-3.0, 5.0), "did not get expected minmax"
|
|
|
|
|
|
###############################################################################
|
|
# Test read and copy of dataset with nan as nodata value (#3576)
|
|
|
|
|
|
def test_stats_nan_3():
|
|
|
|
src_ds = gdal.Open("data/nan32_nodata.tif")
|
|
nodata = src_ds.GetRasterBand(1).GetNoDataValue()
|
|
assert gdaltest.isnan(nodata), "expected nan, got %f" % nodata
|
|
|
|
out_ds = gdaltest.gtiff_drv.CreateCopy("tmp/nan32_nodata.tif", src_ds)
|
|
del out_ds
|
|
|
|
src_ds = None
|
|
|
|
try:
|
|
os.remove("tmp/nan32_nodata.tif.aux.xml")
|
|
except OSError:
|
|
pass
|
|
|
|
ds = gdal.Open("tmp/nan32_nodata.tif")
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
gdaltest.gtiff_drv.Delete("tmp/nan32_nodata.tif")
|
|
assert gdaltest.isnan(nodata), "expected nan, got %f" % nodata
|
|
|
|
|
|
###############################################################################
|
|
# Test reading a VRT with a complex source that define nan as band nodata
|
|
# and complex source nodata (#3576)
|
|
|
|
|
|
def test_stats_nan_4():
|
|
|
|
ds = gdal.Open("data/nan32_nodata.vrt")
|
|
cs = ds.GetRasterBand(1).Checksum()
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
assert cs == 874, "did not get expected checksum"
|
|
|
|
assert gdaltest.isnan(nodata), "expected nan, got %f" % nodata
|
|
|
|
|
|
###############################################################################
|
|
# Test reading a VRT with a complex source that define 0 as band nodata
|
|
# and complex source nodata (nan must be translated to 0 then) (#3576)
|
|
|
|
|
|
def test_stats_nan_5():
|
|
|
|
ds = gdal.Open("data/nan32_nodata_nan_to_zero.vrt")
|
|
cs = ds.GetRasterBand(1).Checksum()
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
assert cs == 978, "did not get expected checksum"
|
|
|
|
assert nodata == 0, "expected nan, got %f" % nodata
|
|
|
|
|
|
###############################################################################
|
|
# Test reading a warped VRT with nan as src nodata and dest nodata (#3576)
|
|
|
|
|
|
def test_stats_nan_6():
|
|
|
|
ds = gdal.Open("data/nan32_nodata_warp.vrt")
|
|
cs = ds.GetRasterBand(1).Checksum()
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
assert cs == 874, "did not get expected checksum"
|
|
|
|
assert gdaltest.isnan(nodata), "expected nan, got %f" % nodata
|
|
|
|
|
|
###############################################################################
|
|
# Test reading a warped VRT with nan as src nodata and 0 as dest nodata (#3576)
|
|
|
|
|
|
def test_stats_nan_7():
|
|
|
|
ds = gdal.Open("data/nan32_nodata_warp_nan_to_zero.vrt")
|
|
cs = ds.GetRasterBand(1).Checksum()
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
assert cs == 978, "did not get expected checksum"
|
|
|
|
assert nodata == 0, "expected nan, got %f" % nodata
|
|
|
|
|
|
###############################################################################
|
|
# Test reading a warped VRT with zero as src nodata and nan as dest nodata (#3576)
|
|
|
|
|
|
def test_stats_nan_8():
|
|
|
|
ds = gdal.Open("data/nan32_nodata_warp_zero_to_nan.vrt")
|
|
cs = ds.GetRasterBand(1).Checksum()
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
assert cs == 874, "did not get expected checksum"
|
|
|
|
assert gdaltest.isnan(nodata), "expected nan, got %f" % nodata
|
|
|
|
|
|
###############################################################################
|
|
# Test statistics computation when nodata = +/- inf
|
|
|
|
|
|
def stats_nodata_inf_progress_cbk(value, string, extra):
|
|
# pylint: disable=unused-argument
|
|
extra[0] = value
|
|
|
|
|
|
def test_stats_nodata_inf():
|
|
|
|
ds = gdal.GetDriverByName("HFA").Create(
|
|
"/vsimem/stats_nodata_inf.img", 3, 1, 1, gdal.GDT_Float32
|
|
)
|
|
ds.GetRasterBand(1).SetNoDataValue(gdaltest.neginf())
|
|
ds.GetRasterBand(1).WriteRaster(
|
|
0, 0, 1, 1, struct.pack("f", gdaltest.neginf()), buf_type=gdal.GDT_Float32
|
|
)
|
|
ds.GetRasterBand(1).WriteRaster(
|
|
1, 0, 1, 1, struct.pack("f", 1), buf_type=gdal.GDT_Float32
|
|
)
|
|
ds.GetRasterBand(1).WriteRaster(
|
|
2, 0, 1, 1, struct.pack("f", -2), buf_type=gdal.GDT_Float32
|
|
)
|
|
|
|
ds.GetRasterBand(1).Checksum()
|
|
user_data = [0]
|
|
stats = ds.GetRasterBand(1).ComputeStatistics(
|
|
False, stats_nodata_inf_progress_cbk, user_data
|
|
)
|
|
assert user_data[0] == 1.0, "did not get expected pct"
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("HFA").Delete("/vsimem/stats_nodata_inf.img")
|
|
|
|
assert stats == [-2.0, 1.0, -0.5, 1.5], "did not get expected stats"
|
|
|
|
|
|
###############################################################################
|
|
# Test deserialization of +inf/-inf values written by Linux and Windows
|
|
|
|
|
|
def stats_nodata_check(filename, expected_nodata):
|
|
ds = gdal.Open(filename)
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
assert nodata == expected_nodata, "did not get expected nodata value"
|
|
|
|
|
|
def test_stats_nodata_neginf_linux():
|
|
return stats_nodata_check("data/stats_nodata_neginf.tif", gdaltest.neginf())
|
|
|
|
|
|
def test_stats_nodata_neginf_msvc():
|
|
return stats_nodata_check("data/stats_nodata_neginf_msvc.tif", gdaltest.neginf())
|
|
|
|
|
|
def test_stats_nodata_posinf_linux():
|
|
return stats_nodata_check("data/stats_nodata_posinf.tif", gdaltest.posinf())
|
|
|
|
|
|
def test_stats_nodata_posinf_msvc():
|
|
return stats_nodata_check("data/stats_nodata_posinf_msvc.tif", gdaltest.posinf())
|
|
|
|
|
|
###############################################################################
|
|
# Test standard deviation computation on huge values
|
|
|
|
|
|
@pytest.mark.require_driver("AAIGRID")
|
|
def test_stats_stddev_huge_values():
|
|
|
|
gdal.FileFromMemBuffer(
|
|
"/vsimem/stats_stddev_huge_values.asc",
|
|
"""ncols 4
|
|
nrows 4
|
|
xllcorner 0
|
|
yllcorner 0
|
|
cellsize 1
|
|
100000000 100000002 100000000 100000002
|
|
100000000 100000002 100000000 100000002
|
|
100000000 100000002 100000000 100000002
|
|
100000000 100000002 100000000 100000002""",
|
|
)
|
|
ds = gdal.Open("/vsimem/stats_stddev_huge_values.asc")
|
|
stats = ds.GetRasterBand(1).ComputeStatistics(0)
|
|
assert stats == [
|
|
100000000.0,
|
|
100000002.0,
|
|
100000001.0,
|
|
1.0,
|
|
], "did not get expected stats"
|
|
ds = None
|
|
gdal.GetDriverByName("AAIGRID").Delete("/vsimem/stats_stddev_huge_values.asc")
|
|
|
|
|
|
###############################################################################
|
|
# Test approximate statistics computation on a square shaped raster whose first column
|
|
# of blocks is nodata only
|
|
|
|
|
|
def test_stats_square_shape():
|
|
|
|
ds = gdal.GetDriverByName("GTiff").Create(
|
|
"/vsimem/stats_square_shape.tif",
|
|
32,
|
|
32,
|
|
options=["TILED=YES", "BLOCKXSIZE=16", "BLOCKYSIZE=16"],
|
|
)
|
|
ds.GetRasterBand(1).SetNoDataValue(0)
|
|
ds.GetRasterBand(1).WriteRaster(
|
|
16, 0, 16, 32, struct.pack("B" * 1, 255), buf_xsize=1, buf_ysize=1
|
|
)
|
|
stats = ds.GetRasterBand(1).ComputeStatistics(True)
|
|
hist = ds.GetRasterBand(1).GetHistogram(approx_ok=1)
|
|
minmax = ds.GetRasterBand(1).ComputeRasterMinMax(1)
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("/vsimem/stats_square_shape.tif")
|
|
|
|
assert stats == [255, 255, 255, 0], "did not get expected stats"
|
|
assert hist[255] == 16 * 16, "did not get expected histogram"
|
|
if minmax != (255, 255):
|
|
print(hist)
|
|
pytest.fail("did not get expected minmax")
|
|
|
|
|
|
###############################################################################
|
|
# Test when nodata = FLT_MIN (#6578)
|
|
|
|
|
|
def test_stats_flt_min():
|
|
|
|
shutil.copyfile("data/flt_min.tif", "tmp/flt_min.tif")
|
|
try:
|
|
os.remove("tmp/flt_min.tif.aux.xml")
|
|
except OSError:
|
|
pass
|
|
|
|
ds = gdal.Open("tmp/flt_min.tif")
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
os.remove("tmp/flt_min.tif.aux.xml")
|
|
|
|
ds = gdal.Open("tmp/flt_min.tif")
|
|
minmax = ds.GetRasterBand(1).ComputeRasterMinMax()
|
|
ds = None
|
|
|
|
os.remove("tmp/flt_min.tif")
|
|
|
|
if nodata != 1.17549435082228751e-38:
|
|
print("%.18g" % nodata)
|
|
pytest.fail("did not get expected nodata")
|
|
|
|
assert stats == [0.0, 1.0, 0.33333333333333337, 0.47140452079103168] or stats == [
|
|
0.0,
|
|
1.0,
|
|
0.33333333333333331,
|
|
0.47140452079103168,
|
|
], "did not get expected stats"
|
|
|
|
assert minmax == (0.0, 1.0), "did not get expected minmax"
|
|
|
|
|
|
###############################################################################
|
|
# Test when nodata = DBL_MIN (#6578)
|
|
|
|
|
|
def test_stats_dbl_min():
|
|
|
|
shutil.copyfile("data/dbl_min.tif", "tmp/dbl_min.tif")
|
|
try:
|
|
os.remove("tmp/dbl_min.tif.aux.xml")
|
|
except OSError:
|
|
pass
|
|
|
|
ds = gdal.Open("tmp/dbl_min.tif")
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
nodata = ds.GetRasterBand(1).GetNoDataValue()
|
|
ds = None
|
|
|
|
os.remove("tmp/dbl_min.tif.aux.xml")
|
|
|
|
ds = gdal.Open("tmp/dbl_min.tif")
|
|
minmax = ds.GetRasterBand(1).ComputeRasterMinMax()
|
|
ds = None
|
|
|
|
os.remove("tmp/dbl_min.tif")
|
|
|
|
if nodata != 2.22507385850720138e-308:
|
|
print("%.18g" % nodata)
|
|
pytest.fail("did not get expected nodata")
|
|
|
|
assert stats == [0.0, 1.0, 0.33333333333333337, 0.47140452079103168] or stats == [
|
|
0.0,
|
|
1.0,
|
|
0.33333333333333331,
|
|
0.47140452079103168,
|
|
], "did not get expected stats"
|
|
|
|
assert minmax == (0.0, 1.0), "did not get expected minmax"
|
|
|
|
|
|
###############################################################################
|
|
# Test stats on a tiled Byte with partial tiles
|
|
|
|
|
|
def test_stats_byte_partial_tiles():
|
|
|
|
ds = gdal.Translate(
|
|
"/vsimem/stats_byte_tiled.tif",
|
|
"../gdrivers/data/small_world.tif",
|
|
creationOptions=["TILED=YES", "BLOCKXSIZE=64", "BLOCKYSIZE=64"],
|
|
)
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("/vsimem/stats_byte_tiled.tif")
|
|
|
|
expected_stats = [0.0, 255.0, 50.22115, 67.119029288849973]
|
|
assert stats == expected_stats, "did not get expected stats"
|
|
|
|
# Same but with nodata set
|
|
ds = gdal.Translate(
|
|
"/vsimem/stats_byte_tiled.tif",
|
|
"../gdrivers/data/small_world.tif",
|
|
creationOptions=["TILED=YES", "BLOCKXSIZE=64", "BLOCKYSIZE=64"],
|
|
)
|
|
ds.GetRasterBand(1).SetNoDataValue(0)
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("/vsimem/stats_byte_tiled.tif")
|
|
|
|
expected_stats = [1.0, 255.0, 50.311081057390084, 67.14541389488096]
|
|
assert stats == pytest.approx(
|
|
expected_stats, rel=1e-10
|
|
), "did not get expected stats"
|
|
|
|
# Same but with nodata set but untiled and with non power of 16 block size
|
|
ds = gdal.Translate(
|
|
"/vsimem/stats_byte_untiled.tif",
|
|
"../gdrivers/data/small_world.tif",
|
|
options="-srcwin 0 0 399 200",
|
|
)
|
|
ds.GetRasterBand(1).SetNoDataValue(0)
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("/vsimem/stats_byte_untiled.tif")
|
|
|
|
expected_stats = [1.0, 255.0, 50.378183963744554, 67.184793517649453]
|
|
assert stats == pytest.approx(
|
|
expected_stats, rel=1e-10
|
|
), "did not get expected stats"
|
|
|
|
ds = gdal.GetDriverByName("GTiff").Create(
|
|
"/vsimem/stats_byte_tiled.tif",
|
|
1000,
|
|
512,
|
|
options=["TILED=YES", "BLOCKXSIZE=512", "BLOCKYSIZE=512"],
|
|
)
|
|
ds.GetRasterBand(1).Fill(255)
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
gdal.Unlink("/vsimem/stats_byte_tiled.tif")
|
|
|
|
expected_stats = [255.0, 255.0, 255.0, 0.0]
|
|
assert (
|
|
max([abs(stats[i] - expected_stats[i]) for i in range(4)]) <= 1e-15
|
|
), "did not get expected stats"
|
|
|
|
# Non optimized code path
|
|
ds = gdal.GetDriverByName("MEM").Create("", 1, 1)
|
|
ds.GetRasterBand(1).WriteRaster(0, 0, 1, 1, struct.pack("B" * 1, 1))
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
expected_stats = [1.0, 1.0, 1.0, 0.0]
|
|
assert (
|
|
max([abs(stats[i] - expected_stats[i]) for i in range(4)]) <= 1e-15
|
|
), "did not get expected stats"
|
|
|
|
ds = gdal.GetDriverByName("MEM").Create("", 3, 5)
|
|
ds.GetRasterBand(1).WriteRaster(0, 0, 3, 1, struct.pack("B" * 3, 20, 30, 50))
|
|
ds.GetRasterBand(1).WriteRaster(0, 1, 3, 1, struct.pack("B" * 3, 60, 10, 5))
|
|
ds.GetRasterBand(1).WriteRaster(0, 2, 3, 1, struct.pack("B" * 3, 10, 20, 0))
|
|
ds.GetRasterBand(1).WriteRaster(0, 3, 3, 1, struct.pack("B" * 3, 10, 20, 255))
|
|
ds.GetRasterBand(1).WriteRaster(0, 4, 3, 1, struct.pack("B" * 3, 10, 20, 10))
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
expected_stats = [0.0, 255.0, 35.333333333333336, 60.785597709398971]
|
|
assert (
|
|
max([abs(stats[i] - expected_stats[i]) for i in range(4)]) <= 1e-15
|
|
), "did not get expected stats"
|
|
|
|
ds = gdal.GetDriverByName("MEM").Create("", 32 + 2, 2)
|
|
ds.GetRasterBand(1).Fill(1)
|
|
ds.GetRasterBand(1).WriteRaster(32, 1, 2, 1, struct.pack("B" * 2, 0, 255))
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
expected_stats = [0.0, 255.0, 4.7205882352941178, 30.576733555893391]
|
|
assert (
|
|
max([abs(stats[i] - expected_stats[i]) for i in range(4)]) <= 1e-15
|
|
), "did not get expected stats"
|
|
|
|
ds = gdal.GetDriverByName("MEM").Create("", 32 + 2, 2)
|
|
ds.GetRasterBand(1).Fill(1)
|
|
ds.GetRasterBand(1).SetNoDataValue(2)
|
|
ds.GetRasterBand(1).WriteRaster(32, 1, 2, 1, struct.pack("B" * 2, 0, 255))
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
expected_stats = [0.0, 255.0, 4.7205882352941178, 30.576733555893391]
|
|
assert (
|
|
max([abs(stats[i] - expected_stats[i]) for i in range(4)]) <= 1e-15
|
|
), "did not get expected stats"
|
|
|
|
|
|
###############################################################################
|
|
# Test stats on uint16
|
|
|
|
|
|
def test_stats_uint16():
|
|
|
|
ds = gdal.Translate(
|
|
"/vsimem/stats_uint16_tiled.tif",
|
|
"../gdrivers/data/small_world.tif",
|
|
outputType=gdal.GDT_UInt16,
|
|
scaleParams=[[0, 255, 0, 65535]],
|
|
creationOptions=["TILED=YES", "BLOCKXSIZE=64", "BLOCKYSIZE=64"],
|
|
)
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("/vsimem/stats_uint16_tiled.tif")
|
|
|
|
expected_stats = [
|
|
0.0,
|
|
65535.0,
|
|
50.22115 * 65535 / 255,
|
|
67.119029288849973 * 65535 / 255,
|
|
]
|
|
assert stats == expected_stats, "did not get expected stats"
|
|
|
|
ds = gdal.Translate(
|
|
"/vsimem/stats_uint16_untiled.tif",
|
|
"../gdrivers/data/small_world.tif",
|
|
options="-srcwin 0 0 399 200 -scale 0 255 0 65535 -ot UInt16",
|
|
)
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("/vsimem/stats_uint16_untiled.tif")
|
|
|
|
expected_stats = [0.0, 65535.0, 12923.9921679198, 17259.703026841547]
|
|
assert stats == expected_stats, "did not get expected stats"
|
|
|
|
# Same but with nodata set but untiled and with non power of 16 block size
|
|
ds = gdal.Translate(
|
|
"/vsimem/stats_uint16_untiled.tif",
|
|
"../gdrivers/data/small_world.tif",
|
|
options="-srcwin 0 0 399 200 -scale 0 255 0 65535 -ot UInt16",
|
|
)
|
|
ds.GetRasterBand(1).SetNoDataValue(0)
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("GTiff").Delete("/vsimem/stats_uint16_untiled.tif")
|
|
|
|
expected_stats = [
|
|
257.0,
|
|
65535.0,
|
|
50.378183963744554 * 65535 / 255,
|
|
67.184793517649453 * 65535 / 255,
|
|
]
|
|
assert stats == pytest.approx(
|
|
expected_stats, rel=1e-10
|
|
), "did not get expected stats"
|
|
|
|
for fill_val in [0, 1, 32767, 32768, 65535]:
|
|
ds = gdal.GetDriverByName("GTiff").Create(
|
|
"/vsimem/stats_uint16_tiled.tif",
|
|
1000,
|
|
512,
|
|
1,
|
|
gdal.GDT_UInt16,
|
|
options=["TILED=YES", "BLOCKXSIZE=512", "BLOCKYSIZE=512"],
|
|
)
|
|
ds.GetRasterBand(1).Fill(fill_val)
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
gdal.Unlink("/vsimem/stats_uint16_tiled.tif")
|
|
|
|
expected_stats = [fill_val, fill_val, fill_val, 0.0]
|
|
if max([abs(stats[i] - expected_stats[i]) for i in range(4)]) > 1e-15:
|
|
print(fill_val)
|
|
pytest.fail("did not get expected stats")
|
|
|
|
# Test remaining pixels after multiple of 32
|
|
ds = gdal.GetDriverByName("MEM").Create("", 32 + 2, 1, 1, gdal.GDT_UInt16)
|
|
ds.GetRasterBand(1).Fill(1)
|
|
ds.GetRasterBand(1).WriteRaster(32, 0, 2, 1, struct.pack("H" * 2, 0, 65535))
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
expected_stats = [0.0, 65535.0, 1928.4411764705883, 11072.48066469611]
|
|
assert (
|
|
max([abs(stats[i] - expected_stats[i]) for i in range(4)]) <= 1e-15
|
|
), "did not get expected stats"
|
|
|
|
# Non optimized code path
|
|
for fill_val in [0, 1, 32767, 32768, 65535]:
|
|
ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1, gdal.GDT_UInt16)
|
|
ds.GetRasterBand(1).WriteRaster(0, 0, 1, 1, struct.pack("H" * 1, fill_val))
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
expected_stats = [fill_val, fill_val, fill_val, 0.0]
|
|
if max([abs(stats[i] - expected_stats[i]) for i in range(4)]) > 1e-15:
|
|
print(fill_val)
|
|
pytest.fail("did not get expected stats")
|
|
|
|
ds = gdal.GetDriverByName("MEM").Create("", 3, 5, 1, gdal.GDT_UInt16)
|
|
ds.GetRasterBand(1).WriteRaster(0, 0, 3, 1, struct.pack("H" * 3, 20, 30, 50))
|
|
ds.GetRasterBand(1).WriteRaster(0, 1, 3, 1, struct.pack("H" * 3, 60, 10, 5))
|
|
ds.GetRasterBand(1).WriteRaster(0, 2, 3, 1, struct.pack("H" * 3, 10, 20, 0))
|
|
ds.GetRasterBand(1).WriteRaster(0, 3, 3, 1, struct.pack("H" * 3, 10, 20, 65535))
|
|
ds.GetRasterBand(1).WriteRaster(0, 4, 3, 1, struct.pack("H" * 3, 10, 20, 10))
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
expected_stats = [0.0, 65535.0, 4387.333333333333, 16342.408927558861]
|
|
assert (
|
|
max([abs(stats[i] - expected_stats[i]) for i in range(4)]) <= 1e-15
|
|
), "did not get expected stats"
|
|
|
|
ds = gdal.GetDriverByName("MEM").Create("", 2, 2, 1, gdal.GDT_UInt16)
|
|
ds.GetRasterBand(1).WriteRaster(0, 0, 2, 1, struct.pack("H" * 2, 0, 65535))
|
|
ds.GetRasterBand(1).WriteRaster(0, 1, 2, 1, struct.pack("H" * 2, 1, 65534))
|
|
stats = ds.GetRasterBand(1).GetStatistics(0, 1)
|
|
ds = None
|
|
|
|
expected_stats = [0.0, 65535.0, 32767.5, 32767.000003814814]
|
|
assert (
|
|
max([abs(stats[i] - expected_stats[i]) for i in range(4)]) <= 1e-15
|
|
), "did not get expected stats"
|
|
|
|
|
|
###############################################################################
|
|
# Test a case where the nodata value is almost the maximum value of float32
|
|
|
|
|
|
def test_stats_nodata_almost_max_float32():
|
|
|
|
gdal.FileFromMemBuffer(
|
|
"/vsimem/float32_almost_nodata_max_float32.tif",
|
|
open("data/float32_almost_nodata_max_float32.tif", "rb").read(),
|
|
)
|
|
|
|
ds = gdal.Open("/vsimem/float32_almost_nodata_max_float32.tif")
|
|
minmax = ds.GetRasterBand(1).ComputeRasterMinMax()
|
|
assert minmax == (0, 0), "did not get expected minmax"
|
|
stats = ds.GetRasterBand(1).ComputeStatistics(False)
|
|
assert stats == [0, 0, 0, 0], "did not get expected stats"
|
|
hist = ds.GetRasterBand(1).GetHistogram(approx_ok=0)
|
|
assert hist[0] == 3, "did not get expected hist"
|
|
ds = None
|
|
|
|
gdal.GetDriverByName("GTiff").Delete(
|
|
"/vsimem/float32_almost_nodata_max_float32.tif"
|
|
)
|
|
|
|
|
|
###############################################################################
|
|
# Test STATISTICS_APPROXIMATE
|
|
|
|
|
|
def test_stats_approx_stats_flag(dt=gdal.GDT_Byte, struct_frmt="B"):
|
|
|
|
ds = gdal.GetDriverByName("MEM").Create("", 2000, 2000, 1, dt)
|
|
ds.GetRasterBand(1).WriteRaster(1000, 1000, 1, 1, struct.pack(struct_frmt * 1, 20))
|
|
|
|
approx_ok = 1
|
|
force = 1
|
|
stats = ds.GetRasterBand(1).GetStatistics(approx_ok, force)
|
|
assert stats == [0.0, 0.0, 0.0, 0.0], "did not get expected stats"
|
|
md = ds.GetRasterBand(1).GetMetadata()
|
|
assert md == {
|
|
"STATISTICS_MEAN": "0",
|
|
"STATISTICS_MAXIMUM": "0",
|
|
"STATISTICS_MINIMUM": "0",
|
|
"STATISTICS_APPROXIMATE": "YES",
|
|
"STATISTICS_STDDEV": "0",
|
|
"STATISTICS_VALID_PERCENT": "100",
|
|
}, "did not get expected metadata"
|
|
|
|
approx_ok = 0
|
|
force = 0
|
|
stats = ds.GetRasterBand(1).GetStatistics(approx_ok, force)
|
|
assert stats == [0.0, 0.0, 0.0, -1.0], "did not get expected stats"
|
|
|
|
approx_ok = 0
|
|
force = 1
|
|
stats = ds.GetRasterBand(1).GetStatistics(approx_ok, force)
|
|
assert stats[1] == 20.0, "did not get expected stats"
|
|
md = ds.GetRasterBand(1).GetMetadata()
|
|
assert "STATISTICS_APPROXIMATE" not in md, "did not get expected metadata"
|
|
|
|
|
|
def test_stats_approx_stats_flag_float():
|
|
return test_stats_approx_stats_flag(dt=gdal.GDT_Float32, struct_frmt="f")
|
|
|
|
|
|
def test_stats_all_nodata():
|
|
|
|
ds = gdal.GetDriverByName("MEM").Create("", 2000, 2000)
|
|
ds.GetRasterBand(1).SetNoDataValue(0)
|
|
|
|
with gdaltest.disable_exceptions(), gdaltest.error_handler():
|
|
minmax = ds.GetRasterBand(1).ComputeRasterMinMax()
|
|
assert math.isnan(minmax[0])
|
|
assert math.isnan(minmax[1])
|
|
|
|
with pytest.raises(Exception):
|
|
ds.GetRasterBand(1).ComputeRasterMinMax()
|
|
|
|
with pytest.raises(Exception):
|
|
ds.GetRasterBand(1).ComputeRasterMinMax(can_return_none=True)
|
|
|
|
with pytest.raises(Exception):
|
|
# can_return_null also accepted for similarity with other methods
|
|
ds.GetRasterBand(1).ComputeRasterMinMax(can_return_null=True)
|
|
|
|
approx_ok = 1
|
|
force = 1
|
|
with pytest.raises(Exception):
|
|
ds.GetRasterBand(1).GetStatistics(approx_ok, force)
|
|
|
|
ds = gdal.GetDriverByName("MEM").Create("", 2000, 2000, 1, gdal.GDT_Float32)
|
|
ds.GetRasterBand(1).SetNoDataValue(0)
|
|
approx_ok = 1
|
|
force = 1
|
|
with pytest.raises(Exception):
|
|
ds.GetRasterBand(1).GetStatistics(approx_ok, force)
|
|
|
|
|
|
def test_stats_float32_with_nodata_slightly_above_float_max():
|
|
|
|
ds = gdal.Open("data/float32_with_nodata_slightly_above_float_max.tif")
|
|
my_min, my_max = ds.GetRasterBand(1).ComputeRasterMinMax()
|
|
assert (my_min, my_max) == (
|
|
-1.0989999771118164,
|
|
0.703338623046875,
|
|
), "did not get expected stats"
|
|
|
|
|
|
def test_stats_clear():
|
|
|
|
filename = "/vsimem/out.tif"
|
|
gdal.Translate(filename, "data/byte.tif")
|
|
ds = gdal.Open(filename)
|
|
assert ds.GetRasterBand(1).GetStatistics(False, False) == [0, 0, 0, -1]
|
|
assert ds.GetRasterBand(1).ComputeStatistics(False) != [0, 0, 0, -1]
|
|
|
|
ds = gdal.Open(filename)
|
|
assert ds.GetRasterBand(1).GetStatistics(False, False) != [0, 0, 0, -1]
|
|
ds.ClearStatistics()
|
|
|
|
ds = gdal.Open(filename)
|
|
assert ds.GetRasterBand(1).GetStatistics(False, False) == [0, 0, 0, -1]
|
|
|
|
gdal.GetDriverByName("GTiff").Delete(filename)
|
|
|
|
|
|
@pytest.mark.parametrize(
|
|
"datatype,minval,maxval",
|
|
[
|
|
(gdal.GDT_Byte, 1, 254),
|
|
(gdal.GDT_Byte, -127, 127),
|
|
(gdal.GDT_UInt16, 1, 65535),
|
|
(gdal.GDT_Int16, -32767, 32766),
|
|
(gdal.GDT_UInt32, 1, (1 << 32) - 2),
|
|
(gdal.GDT_Int32, -(1 << 31) + 1, (1 << 31) - 2),
|
|
(gdal.GDT_UInt64, 1, (1 << 53) - 2),
|
|
(gdal.GDT_Int64, -(1 << 53) + 2, (1 << 53) - 2),
|
|
(
|
|
gdal.GDT_Float32,
|
|
-struct.unpack("f", struct.pack("f", 1e20))[0],
|
|
struct.unpack("f", struct.pack("f", 1e20))[0],
|
|
),
|
|
(gdal.GDT_Float64, -1e100, 1e100),
|
|
(gdal.GDT_CInt16, -32767, 32766),
|
|
(gdal.GDT_CInt32, -(1 << 31) + 1, (1 << 31) - 2),
|
|
(
|
|
gdal.GDT_CFloat32,
|
|
-struct.unpack("f", struct.pack("f", 1e20))[0],
|
|
struct.unpack("f", struct.pack("f", 1e20))[0],
|
|
),
|
|
(gdal.GDT_CFloat64, -1e100, 1e100),
|
|
],
|
|
)
|
|
@pytest.mark.parametrize("nodata", [None, 0, 1])
|
|
def test_stats_computeminmax(datatype, minval, maxval, nodata):
|
|
ds = gdal.GetDriverByName("MEM").Create("", 64, 1, 1, datatype)
|
|
minval_mod = minval
|
|
expected_minval = minval
|
|
if datatype == gdal.GDT_Byte and minval < 0:
|
|
ds.GetRasterBand(1).SetMetadataItem(
|
|
"PIXELTYPE", "SIGNEDBYTE", "IMAGE_STRUCTURE"
|
|
)
|
|
minval_mod = 256 + minval
|
|
if nodata:
|
|
ds.GetRasterBand(1).SetNoDataValue(nodata)
|
|
if nodata == 1 and minval == 1:
|
|
expected_minval = maxval
|
|
ds.GetRasterBand(1).WriteRaster(
|
|
0,
|
|
0,
|
|
64,
|
|
1,
|
|
struct.pack("d" * 2, minval_mod, maxval),
|
|
buf_type=gdal.GDT_Float64,
|
|
buf_xsize=2,
|
|
buf_ysize=1,
|
|
)
|
|
assert ds.GetRasterBand(1).ComputeRasterMinMax(0) == (expected_minval, maxval)
|
|
|
|
|
|
###############################################################################
|
|
# Test statistics on band with mask band
|
|
|
|
|
|
def test_stats_mask_band():
|
|
|
|
src_ds = gdal.GetDriverByName("MEM").Create("", 3, 1, 1, gdal.GDT_Int16)
|
|
src_ds.WriteRaster(0, 0, 3, 1, struct.pack("h" * 3, 1, 2, 3))
|
|
src_ds.CreateMaskBand(gdal.GMF_PER_DATASET)
|
|
mask_band = src_ds.GetRasterBand(1).GetMaskBand()
|
|
mask_band.WriteRaster(0, 0, 3, 1, struct.pack("B" * 3, 0, 255, 255))
|
|
assert src_ds.GetRasterBand(1).ComputeRasterMinMax(False) == (2, 3)
|
|
assert src_ds.GetRasterBand(1).ComputeStatistics(False) == [2, 3, 2.5, 0.5]
|
|
assert src_ds.GetRasterBand(1).GetHistogram(False) == [0, 0, 1, 1] + ([0] * 252)
|