TransWikia.com

Black and white to rgb

Geographic Information Systems Asked on June 4, 2021

I have a GeoTIFF grayscale raster. When trying to paint a raster in RGB, the code does not take into account the nodata value.

import sys
from osgeo import gdal, gdalnumeric
import sys
from xml.dom import minidom
import numpy as np
from os.path import join
import logging
from baseimage import BaseImage
import glob, os

dirName = sys.argv[1]

gdal.UseExceptions()

os.chdir("/var/www/html/algorithm/tiles/base_tiles/cut_38ULV_prod")

class ColorMap(BaseImage):
    """
    Class parse color map qml file and apply color map to GDAL data set.
    """
    MULTIPL = 1000
    STEP = 1
    MINVAL = 0
    MAXVAL = 255
      
    def __init__(self, color_map_path, nodata):
        
        self.color_map_path = color_map_path
        self.nodata = nodata
        self.cm = {}
        self.min = None
        self.max = None
        
        
    def _rgb(self, triplet):
        _NUMERALS = '0123456789abcdefABCDEF'
        _HEXDEC = {v: int(v, 16) for v in (x+y for x in _NUMERALS for y in _NUMERALS)}
    
        return [_HEXDEC[triplet[0:2]], _HEXDEC[triplet[2:4]], _HEXDEC[triplet[4:6]]]
        
    
    def _rescale_formula(self, initial_value, _initial_max, _initial_min, _rescale_max, _rescale_min, _reverse):      
        _initial_interval = _initial_max - _initial_min
        _rescale_interval = _rescale_max - _rescale_min
                
        if _reverse: 
            return round(_rescale_max - (float((initial_value - _initial_min) * _rescale_interval)/float(_initial_interval))) 
        else:
            return round((float((initial_value - _initial_min) * _rescale_interval)/float(_initial_interval)) + _rescale_min) 
        
    def _rescale_band(self, _initial_max, _initial_min, _rescale_end, _rescale_start, initial_value): 
        _rescale_nodata = self.MINVAL
        _initial_nodata = self.nodata
        
        if initial_value == _initial_nodata:            
            return _rescale_nodata
        elif _rescale_end == _rescale_start: 
            return _rescale_end
        else:
            _rescale_min = min(_rescale_end, _rescale_start)
            _rescale_max = max(_rescale_end, _rescale_start)
            _reverse = False if _rescale_min == _rescale_start else True
            return self._rescale_formula(initial_value, _initial_max, _initial_min, _rescale_max, _rescale_min, _reverse)
    
    def _get_rgb(self, value, band):        
        if value < self.min: 
            return self.cm[self.min][band]
        elif value > self.max: 
            return self.cm[self.max][band]           
        return self.cm[value][band] 
           
        
    def parse_color_map(self):
        
        colormap = []
                
        dom = minidom.parse(self.color_map_path)
        colorrampshader = dom.getElementsByTagName("colorrampshader")[0]
        
        for item in colorrampshader.getElementsByTagName("item"):    
            colormap.append({'mult_value': round(float(item.getAttribute('value')) * self.MULTIPL) 
                            ,'rgb_color': self._rgb( item.getAttribute('color')[1:] ) })        
                
        for i in range(0, len(colormap)-1):       
            self.cm[colormap[i]['mult_value']] = colormap[i]['rgb_color']
                                                  
            dif = abs(colormap[i+1]['mult_value'] - colormap[i]['mult_value'])
            if dif > self.STEP:
                for j in range(colormap[i]['mult_value'] + self.STEP, colormap[i+1]['mult_value'], self.STEP):            
                    self.cm[j] = [self._rescale_band(colormap[i+1]['mult_value'], colormap[i]['mult_value'], colormap[i+1]['rgb_color'][0], colormap[i]['rgb_color'][0], j)
                                 ,self._rescale_band(colormap[i+1]['mult_value'], colormap[i]['mult_value'], colormap[i+1]['rgb_color'][1], colormap[i]['rgb_color'][1], j)
                                 ,self._rescale_band(colormap[i+1]['mult_value'], colormap[i]['mult_value'], colormap[i+1]['rgb_color'][2], colormap[i]['rgb_color'][2], j)]
         
        self.cm[colormap[i+1]['mult_value']] = [colormap[i+1]['rgb_color'][0], colormap[i+1]['rgb_color'][1], colormap[i+1]['rgb_color'][2]] 
        
        self.min = min(self.cm.keys())
        self.max = max(self.cm.keys())
        
        self.cm[self.nodata * self.MULTIPL] = [self.MINVAL, self.MINVAL, self.MINVAL]
        
        
    def create_rgba(self, src_path, dst_path):
    
        src_ds = self._open_file(src_path, gdal.GA_ReadOnly)        
        
        # dst_path = join(dst_dir, 'ndvi_img.tif')
         
        try:
            # Prepare output file and dataset.
            dst_ds = self._create_output_ds(src_ds, dst_path, self.FORMAT_GTIFF, 4, gdal.GDT_Byte)
                                   
            for line in range(0, src_ds.RasterYSize): 
                line_array = gdalnumeric.BandReadAsArray(src_ds.GetRasterBand(1)                                                      
                                                        ,yoff=line
                                                        ,win_xsize=src_ds.RasterXSize 
                                                        ,win_ysize=1)

                band_array_unt8 = (np.around(line_array * self.MULTIPL)).astype(int)              
                to_rgb = np.vectorize(self._get_rgb, otypes=[np.uint8])        
                red = to_rgb(band_array_unt8, 0)
                green = to_rgb(band_array_unt8, 1)           
                blue = to_rgb(band_array_unt8, 2)
               
                alpha = np.where(line_array == self.nodata, self.MINVAL, self.MAXVAL)         
               
                # Write bands to the output file.                
                gdalnumeric.BandWriteArray(dst_ds.GetRasterBand(1), red, yoff=line)
                gdalnumeric.BandWriteArray(dst_ds.GetRasterBand(2), green, yoff=line)
                gdalnumeric.BandWriteArray(dst_ds.GetRasterBand(3), blue, yoff=line)
                gdalnumeric.BandWriteArray(dst_ds.GetRasterBand(4), alpha, yoff=line)
         
        except RuntimeError as e:
            print('Unable to create rgda NDVI: "{}"'.format(e))
            sys.exit(-1) 
        finally:
            src_ds = None
            dst_ds = None 
            line_array = None
            band_array_unt8 = None
            red = None
            green = None
            blue = None
            alpha = None
    
        return dst_path


if __name__ == '__main__':
    
    print(logging.DEBUG)
    
    cm = ColorMap(color_map_path='/var/www/html/algorithm/ndvi_style_prod_areas.qml', nodata=0)
    cm.parse_color_map()
    for file in glob.glob("*.tif"):
        fileName = file
        newFileName = '_colored_'+fileName
        print(dirName)
        print(fileName)
        cm.create_rgba(src_path="/var/www/html/algorithm/tiles/base_tiles/"+dirName+'/'+fileName, dst_path='/var/www/html/algorithm/tiles/base_tiles/'+dirName+'/'+newFileName)

    print('coloring Done!')

result

and in the place where the alpha channel should be, a purple fill. Please tell me how to get rid of the fill and achieve an alpha channel.

One Answer

The problem has been resolved. The point was that when creating the original Geotiff, I designated nodata = 0, which was an error, since I work with NDVI. Changing the nodata value in the source file to any other solves the problem.

Answered by Константин on June 4, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP