nsidc / qgreenland

Source code for generating the QGreenland package hosted at https://qgreenland.org/
https://qgreenland.readthedocs.io
Other
36 stars 9 forks source link

Fix geoid and gravity anomaly layers' nodata value in overviews #712

Closed trey-stafford closed 11 months ago

trey-stafford commented 11 months ago

The given nodata value is not set correctly in the overviews. Use -9999 instead, which does get set correctly.

In [2]: from osgeo import gdal

In [3]: ds = gdal.Open('/share/appdata/qgreenland/working-storage/wip-layers/geoid/04-command-build_overviews/ggeoid16.tif')

In [4]: band = ds.GetRasterBand(1);

In [5]: band
Out[5]: <osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x7f0dcf56c0c0> >

In [8]: band.GetOverview(2)
Out[8]: <osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x7f0dcea1a310> >

In [9]: band.GetOverview(2).ReadAsArray()
Out[9]:
array([[1.70141001e+38, 1.70141001e+38, 1.70141001e+38, ...,
        1.70141001e+38, 1.70141001e+38, 1.70141001e+38],
       [1.70141001e+38, 1.70141001e+38, 1.70141001e+38, ...,
        1.70141001e+38, 1.70141001e+38, 1.70141001e+38],
       [1.70141001e+38, 1.70141001e+38, 1.70141001e+38, ...,
        1.70141001e+38, 1.70141001e+38, 1.70141001e+38],
       ...,
       [1.70141001e+38, 1.70141001e+38, 1.70141001e+38, ...,
        1.70141001e+38, 1.70141001e+38, 1.70141001e+38],
       [1.70141001e+38, 1.70141001e+38, 1.70141001e+38, ...,
        1.70141001e+38, 1.70141001e+38, 1.70141001e+38],
       [1.70141001e+38, 1.70141001e+38, 1.70141001e+38, ...,
        1.70141001e+38, 1.70141001e+38, 1.70141001e+38]])

In [20]: band.GetOverview(2).ReadAsArray().dtype
Out[20]: dtype('float64')

In [21]: band.ReadAsArray().dtype
Out[21]: dtype('float64')

In [22]: band.GetOverview(2).ReadAsArray()[0, 0]
Out[22]: 1.701410009187828e+38

In [23]: band.ReadAsArray()[0, 0]
Out[23]: 1.70141e+38

In [24]: overview_nodata = band.GetOverview(2).ReadAsArray()[0, 0]

In [25]: data_nodata = band.ReadAsArray()[0, 0]

In [26]: overview_nodata == data_nodata
Out[26]: False

In [27]: overview_nodata < data_nodata
Out[27]: False

In [28]: overview_nodata > data_nodata
Out[28]: True

In [29]: overview_nodata
Out[29]: 1.701410009187828e+38

In [30]: int(overview_nodata)
Out[30]: 170141000918782798866653488190622531584

In [31]: int(data_nodata)
Out[31]: 170141000000000007096036300471486382080

Description

REPLACE ME WITH A PULL REQUEST DESCRIPTION.

Checklist

If an item on this list is done or not needed, simply check it with [x].

MattF-NSIDC commented 11 months ago

We think this may be a GDAL bug!

https://github.com/OSGeo/gdal/issues/8187