Read the documentation of VRT https://gdal.org/drivers/raster/vrt.html and add color table. Step by step example:
Create test image with gdal_create https://gdal.org/programs/gdal_create.html#gdal-create. Utility is included in GDAL version 3.2 and higher.
gdal_create -of GTiff -outsize 10 10 -bands 1 -burn 2 -ot Float32 float_gray.tif
Create a VRT from ...
It might be worth avoiding use of np.NaN altogether. NaN literally means "not a number", and it cannot be converted to an integer.
There are two ways of doing this, depending on the nature of the data, and what the negative numbers mean in that data (it is the negative values that the script is attempting to convert to np.Nan).
To combine each of the TIFF files into one multiband TIFF file, you should use gdalbuildvrt.
The command would be:
gdalbuildvrt -separate -resolution highest -r cubic combined-raster.vrt ...B01 ...B02 ...B03 ...B04 ...B05 ...B06 ...B07 ...B08 ...B8A ...B09 ...B10 ...B11 ...B12
This would provide a single VRT-raster, which can then either be used directly, ...
Thanks to the extensive help of @mdsumner, this solution works by using the polymer R package to break polygons down into a mesh of triangles, calculating the overlapping triangle segments, and then reassembling the triangles into polygons.
This approach is slower but it seems to generally be robust to these non-noded intersection errors which here stem from ...