3

I have an interpolated raster, which is a masked array, that shows where the data is valid. I want to save it as masked TIFF with rasterio.

To test, here's the JSON file, and the code:

import rasterio, rasterio.io, rasterio.transform
import numpy as np
import json

with open('data.json') as f:
    d, m = json.load(f)

zi = np.ma.masked_array(d, m)
res = 100000
minx = 2282862
miny = 5309393
dtype = 'uint8'
tr = rasterio.transform.from_origin(minx - res / 2, miny - res / 2, res, -res)

raster_params = {
    'driver': 'GTiff',
    'width': zi.shape[1],
    'height': zi.shape[0],
    'count': 1,
    'dtype': dtype,
    'compress': 'LZMA',
    'crs': {'init': 'epsg:3857'},
    'transform': tr
}

with rasterio.Env(GDAL_TIFF_INTERNAL_MASK=True):
    with rasterio.open('raster.tif', 'w', **raster_params) as dataset:
        dataset.write(np.array([zi.astype(dtype)]))
        dataset.write_mask(zi.mask)

As a result, there's one band and it has zeroes instead of nan and no mask.

Output_result

I tried using two layers and writing this way:

...
    'count': 2,
...
        dataset.write(np.array([zi.astype(dtype), zi.mask.astype(dtype)]))

It creates a band that I can choose in QGIS. But it just hides the layer completely. When Band 2 has '1' for data and '0' for no-data, it looks like this:

Output_result

If I do the other way around ('1' for no-data, '0' for data), the result is the same.

Am I writing it correctly to GeoTIFF?

How to visualize it properly in QGIS?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.