1

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

Nikko
  • 581
  • 4
  • 16
  • 1
    Check with gdalinfo if the original data happens to be scaled. – user30184 Jan 29 '23 at 14:31
  • What does it mean to be scaled in this manner? In your answer below, the Scale and scale_factor are the same. – Nikko Jan 30 '23 at 08:05
  • 1
    Scale and scalefactor mean that if the raster pixel has a value 20 it must be multiplied by 0.05 for getting the intensity of radiation as watts per square meter that seems to be what the raster data is presenting. Scale is a standard TIFF metadata and scale_factor is stored as NETCDF metadata into another place of the file. In the conversion from NetCDF into TIFF the scale_factor has been converted into Scale that is the corresponding tiff tag but the original metadata scale_factor is also kept. – user30184 Jan 30 '23 at 08:32
  • Thank you. My previous question would consider the scale and offset: https://gis.stackexchange.com/questions/392438/how-to-mask-and-scale-raster-using-rasterio in the output file. – Nikko Jan 30 '23 at 09:06

1 Answers1

1

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
user30184
  • 65,331
  • 4
  • 65
  • 118