Rasterio: Data get promoted to float64 when using WarpedVRT

[vincentsarago]

This was first raised over https://github.com/cogeotiff/rio-cogeo/issues/85
After some digging, I couldn't find the bug myself. I believe the problem comes from our WarpedVRT implementation but when reading the code I cannot find what's going on.

The behaviour only happens with the file shared by the user https://esip-pangeo-uswest2.s3-us-west-2.amazonaws.com/sciencebase/tiny.tif and I couldn't reproduce it with other files.

Expected behavior and actual behavior.

keep original datatype

Steps to reproduce the problem.

File: https://esip-pangeo-uswest2.s3-us-west-2.amazonaws.com/sciencebase/tiny.tif

import rasterio
from rasterio.vrt import WarpedVRT

src_path = "tiny.tif"

with rasterio.open(src_path) as src_dst:
    print(src_dst.meta)
    with WarpedVRT(src_dst) as vrt:
        print(vrt.meta)

> {'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028230607370965e+38, 'width': 512, 'height': 512, 'count': 1, 'crs': CRS.from_epsg(26911), 'transform': Affine(1.0, 0.0, 179523.99999999822,
       0.0, -1.0, 3824832.0)}
> {'driver': 'VRT', 'dtype': 'float64', 'nodata': -3.4028230607370965e+38, 'width': 512, 'height': 512, 'count': 1, 'crs': CRS.from_epsg(26911), 'transform': Affine(1.0, 0.0, 179523.99999999822,
       0.0, -1.0, 3824832.0)}

Operating system

Mac and Ubuntu

Rasterio version and provenance

1.24 and master

cc @rsignell-usgs

[sgillies]

I'm looking into whether we should be setting the working data type when we call GDALCreateWarpedVRT.

[sgillies]

Interestingly, I can't reproduce the issue on my computer with GDAL 3.0.1. PR #1767 is checking with older versions of GDAL.

[vincentsarago]

I'll check with GDAL 3.0 too

[sgillies]

I can't reproduce this on any version of GDAL @vincentsarago . Can you check that my test is correct? https://github.com/mapbox/rasterio/pull/1767/files#diff-e2937f93425979eec69dffba64475224R397

[vincentsarago]

@sgillies I could just reproduce the bug using the file I shared

[sgillies]

@vincentsarago I looked in GDAL and there's a GDALWarpResolveWorkingDataType function that calls GDALDataTypeUnion or GDALDataTypeUnionWithValue. Those functions determine based on the nodata value that float64 is the right data type for the warped VRT. Without a nodata value, you would get "float32". Can you eliminate the nodata value in your case? Or change it to nan? (though I'm not sure what would happen ther) I'm not sure I want to work around GDAL in this situation.

[vincentsarago]

Interesting

Can you eliminate the nodata value in your case?

I can't or at least I don't want to introduce tests for this specific nodata value.

I looked in GDAL and there's a GDALWarpResolveWorkingDataType function that calls GDALDataTypeUnion or GDALDataTypeUnionWithValue

so in GDALWarpResolveWorkingDataType doc it says

If the working data type is unknown,

The point is the dataset data type is known (from src_dst), so I don't get why GDALCreateWarpedVRT change it

[rouault]

If the working data type is unknown,

This is what is specified with -wt on gdalwarp. In GDALCreateWarpedVRT(), the working data type is guessed from the source and target data types, as well as the source and nodata values. It is likely that a nodata value cannot be exactly represented on a Float32, so a change to Float64 is done. This is a big arguable as a behaviour

[sgillies]

@rouault thanks!

@vincentsarago I'm not going to go against GDAL, but I will add a keyword argument to the WarpedVRT constructor so that rio-cogeo can guarantee data type preservation.

[vincentsarago]

thanks @sgillies

I think there is still something here, which make the nodata value to be considered as float64 instead of a float32

with rasterio.open(src_path) as src_dst:
    ...:     print(src_dst.nodata)
    ...:     print(type(src_dst.nodata))
    ...:
-3.4028230607370965e+38
<class 'float'>

But this might be at data creation ...

[vincentsarago]

ah

In [22]: numpy.array([-3.4028230607370965e+38]).dtype
Out[22]: dtype('float64')

~so the nodata value set in the file is considered as a float64 by numpy! This doesn't really make sense to me how we can have a float32 file with a float64 nodata.~

Edit: This starts to make sense now, the nodata value in https://esip-pangeo-uswest2.s3-us-west-2.amazonaws.com/sciencebase/tiny.tif is recognized as a float64 by gdal/numpy but is also a valid float32 value.

In [31]: numpy.float64(-3.4028230607370965e+38) == numpy.float32(-3.4028230607370965e+38)
Out[31]: True

[rouault]

This doesn't really make sense to me how we can have a float32 file with a float64 nodata.

In the GDAL data model, nodata is exposed as a double whatever the raster type is. In GTiff, the nodata value is stored as an ASCII string, so it is well possible to have a nodata value whose type is not compatible with the raster data type. And there might be also rounding issues due to serialization/deserialization to ASCII. This is a recurring source of headaches

[sgillies]

@vincentsarago as Even said, -3.4028230607370965e+38 can't be represented in float32 even though it is mathematically within the range of min and max values for float32. It's a bad nodata value for float32.

-3.4028235e+38 would be better, maybe. Or maybe passing through string serialization would mess it up as well.

tt = numpy.dtype("float32")
>>> numpy.array([numpy.finfo(tt).min])
array([-3.4028235e+38], dtype=float32)

[vincentsarago]

-3.4028235e+38 would be better, maybe.

it won't change the result
```
In [34]: with rasterio.open(src_path) as src_dst:
...: print(src_dst.meta)
...: with WarpedVRT(src_dst, src_crs="EPSG:4326", nodata=-3.4028235e+38) as vrt:
...: print(vrt.meta)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': None, 'width': 3, 'height': 2, 'count': 1, 'crs': None, 'transform': Affine(100.0, 0.0, 0.0,
0.0, 100.0, 0.0)}
{'driver': 'VRT', 'dtype': 'float64', 'nodata': -3.4028235e+38, 'width': 3, 'height': 2, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(100.0, 0.0, 0.0,
0.0, -100.0, 200.0)}

and more funny thing 

**nodata=3.4** --> NOK

In [37]: with rasterio.open(src_path) as src_dst:
...: print(src_dst.meta)
...: with WarpedVRT(src_dst, src_crs="EPSG:4326", nodata=3.4) as vrt:
...: print(vrt.meta)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': None, 'width': 3, 'height': 2, 'count': 1, 'crs': None, 'transform': Affine(100.0, 0.0, 0.0,
0.0, 100.0, 0.0)}
{'driver': 'VRT', 'dtype': 'float64', 'nodata': 3.4, 'width': 3, 'height': 2, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(100.0, 0.0, 0.0,
0.0, -100.0, 200.0)}

**nodata=3** --> OK

In [39]: with rasterio.open(src_path) as src_dst:
...: print(src_dst.meta)
...: with WarpedVRT(src_dst, src_crs="EPSG:4326", nodata=3.) as vrt:
...: print(vrt.meta)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': None, 'width': 3, 'height': 2, 'count': 1, 'crs': None, 'transform': Affine(100.0, 0.0, 0.0,
0.0, 100.0, 0.0)}
{'driver': 'VRT', 'dtype': 'float32', 'nodata': 3.0, 'width': 3, 'height': 2, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(100.0, 0.0, 0.0,
0.0, -100.0, 200.0)}
```