[GIS] From RGB raster to grayscale

gdalqgisrasterrgb

I have been searching for a way to rasterize shp files (Unexpected output when using gdal_rasterize) and it worked with a little help. Now, my code is first rasterizing the largest sip file and than using those dimensions for further layers. For each layer I am inputting a value (COST) which is used as value for each pixel but it is using a RGB format but I actually only need the files to be in grayscale. The RGB is giving a problem when I enter values higher than 255 => than it subtracts n*255 until the value is in the range of 0-255.

Can anybody help me with changing this code to the one that just uses grayscale and where the cost values will be retained like I want it to?

import os
from osgeo import gdal, ogr
import config3

RASTERIZE_COLOR_FIELD = "__color__"
xmi = 0
xma = 0
ymi = 0
yma = 0
syspath = str()
newmap = str()
cost = 0
outputfilep = str()
source_srs = 0

def one():
    while True:
        try:
            global syspath
            syspath = str(input("Systempath: "))
            break
        except (NameError, SyntaxError, ValueError, TypeError):
            print ""
            print("This is not a string. Please do not forget to type your input between quotation marks.")
            print ""
        except (UnicodeEncodeError):
            print("Please do not use special characters in the systempath.")
            print ""

def three():
    while True:
        try:
            global cost
            cost = int(input("Cost per squared meter for this type of land-use: "))
            break
        except (NameError, ValueError, SyntaxError, UnicodeEncodeError):
            print ""
            print("This is not a numerical value, please re-enter.")
            print ""

def rasterize(pixel_size=1):

    one()

    (shapefilefilepath, shapefilename) = os.path.split(syspath)
    (shapefileshortname, extension) = os.path.splitext(shapefilename) 

    outmap = shapefilefilepath + "/"

    three()  

    orig_data_source = ogr.Open(syspath)
    source_ds = ogr.GetDriverByName("Memory").CopyDataSource(orig_data_source, "")
    source_layer = source_ds.GetLayer(0)
    global source_srs
    source_srs = source_layer.GetSpatialRef()
    x_min, x_max, y_min, y_max = source_layer.GetExtent()
    global xmi
    xmi = x_min
    global ymi
    ymi = y_min
    global xma
    xma = x_max
    global yma
    yma = y_max
    field_def = ogr.FieldDefn(RASTERIZE_COLOR_FIELD, ogr.OFTReal)
    source_layer.CreateField(field_def)
    source_layer_def = source_layer.GetLayerDefn()
    field_index = source_layer_def.GetFieldIndex(RASTERIZE_COLOR_FIELD)
    for feature in source_layer:
        feature.SetField(field_index, cost)
        source_layer.SetFeature(feature)
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)
    global outputfilep
    outputfilep = outmap + shapefileshortname + ".tif"

    target_ds = gdal.GetDriverByName('GTiff').Create(outputfilep, x_res,
            y_res, 3, gdal.GDT_Byte)
    target_ds.SetGeoTransform((
            x_min, pixel_size, 0,
            y_max, 0, -pixel_size,
        ))
    if source_srs:
        target_ds.SetProjection(source_srs.ExportToWkt())
    else:
        target_ds.SetProjection('LOCAL_CS["arbitrary"]')
    err = gdal.RasterizeLayer(target_ds, (1,2,3), source_layer,
            burn_values=(0,0,0),
            options=["ATTRIBUTE=%s" % RASTERIZE_COLOR_FIELD])

def rasterize2(pixel_size=1):
    continue1 = 0

    while continue1 != "STOP":

        while True:
            try:
                continue1 = str(input("Systempath: "))
                break
            except (NameError, SyntaxError, ValueError, TypeError):
                print ""
                print("This is not a string. Please do not forget to type your input between quotation marks.")
                print ""
            except (UnicodeEncodeError):
                print("Please do not use special characters in the systempath.")
                print ""

        while True:
            try:
                cost = int(input("Cost per squared meter for this type of land-use: "))
                break
            except (NameError, ValueError, SyntaxError, UnicodeEncodeError):
                print ""
                print("This is not a numerical value, please re-enter.")
                print ""

        (shapefilefilepath, shapefilename) = os.path.split(continue1)
        (shapefileshortname, extension) = os.path.splitext(shapefilename)
        outmap = shapefilefilepath + "/"


        outputfilep = outmap + shapefileshortname + ".tif"

        orig_data_source = ogr.Open(continue1)
        source_ds = ogr.GetDriverByName("Memory").CopyDataSource(orig_data_source, "")
        source_layer = source_ds.GetLayer(0)
        field_def = ogr.FieldDefn(RASTERIZE_COLOR_FIELD, ogr.OFTReal)
        source_layer.CreateField(field_def)
        source_layer_def = source_layer.GetLayerDefn()
        field_index = source_layer_def.GetFieldIndex(RASTERIZE_COLOR_FIELD)
        # Generate random values for the color field (it's here that the value
        # of the attribute should be used, but you get the idea)
        for feature in source_layer:
            feature.SetField(field_index, cost) # DIT IS AANGEPAST
            source_layer.SetFeature(feature)
        x_res = int((xma - xmi) / pixel_size)
        y_res = int((yma - ymi) / pixel_size)
        target_ds = gdal.GetDriverByName('GTiff').Create(outputfilep, x_res,
                y_res, 3, gdal.GDT_Byte)
        target_ds.SetGeoTransform((
                xmi, pixel_size, 0,
                yma, 0, -pixel_size,
            ))
        if source_srs:
            target_ds.SetProjection(source_srs.ExportToWkt())
        else:
            target_ds.SetProjection('LOCAL_CS["arbitrary"]')
        err = gdal.RasterizeLayer(target_ds, (1,2,3), source_layer, burn_values=(0,0,0), options=["ATTRIBUTE=%s" % RASTERIZE_COLOR_FIELD])

        continue1 = input("Press enter if you want to continue inputting shapefiles. Otherwise type STOP: ")

    if err != 0:
        raise Exception("error rasterizing layer: %s" % err)           

if __name__ == '__main__':
    rasterize()
    rasterize2()

Best Answer

If you want to obtain greyscale from a RGB composite, you can perform a HIS transform (see e.g. i.rgb.his) and then take the resulting intensity channel.

Related Question