Rasterio – Masking Changes Values in Output File

maskingrasterrasterio

I am masking a raster file using rasterio. I am doing this using the mask function. The result of the masking changes the value of the pixels when loaded to GIS.

import fiona
import rasterio
import rasterio.mask

with fiona.open("test_polygon.json", "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

with rasterio.open("exampe.tif") as src:
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta
    
out_meta.update({"driver": "GTiff",
                 "dtype": rasterio.float32,
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

with rasterio.open("swr_masked.tif", "w", **out_meta) as dest:
    dest.write(out_image)

The example files are located HERE.

When loaded to a GIS application like QGIS, the value of the masked data changes. It is like it was changed to another data type like int16 or something. I specified the dtype in writing to be the same with the input float32. But still have the same result.

Result in QGIS:
enter image description here

Best Answer

The original image is scaled:

gdalinfo example.tif
...
Band 1 Block=192x21 Type=Int16, ColorInterp=Gray
  NoData Value=-32768
  Unit Type: W/m^2
  Offset: 0,   Scale:0.0500000007450581
  Metadata:
    add_offset=0
    long_name=Shortwave radiation
    missing_value=-32768
    NETCDF_VARNAME=SWR
    scale_factor=0.050000001
    units=W/m^2
    valid_max=26000
    valid_min=0

QGIS knows how to apply scaling and it divides the raw pixel values by 20 but your rasterio code does not preserve the metadata about scaling.

I do not know how to handle the scaling with rasterio but GDAL seems to do it automatically

gdalwarp -ot float32 -cutline test_polygon.json -crop_to_cutline example.tif out1.tif

enter image description here

gdalinfo out1.tif
...
Band 1 Block=2x1 Type=Float32, ColorInterp=Gray
  NoData Value=-32768
  Unit Type: W/m^2
  Offset: 0,   Scale:0.0500000007450581
  Metadata:
    add_offset=0
    long_name=Shortwave radiation
    missing_value=-32768
    NETCDF_VARNAME=SWR
    scale_factor=0.050000001
    units=W/m^2
    valid_max=26000
    valid_min=0
Related Question