2

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()
  • I don't understand what you're trying to do in Python. – gene May 11 '14 at 15:33
  • What range of values will your "COST" have? – AndreJ May 11 '14 at 15:39
  • The shp files are being rasterized automatically and are getting a value 'cost' for each pixel. These values are representing costs for digging a hole there and will go up to 10000. – Peter Vanvoorden May 11 '14 at 18:16
0

You are not restricted to use RGB raster files with 3 bands ranging from 0 to 255.

The GDAL Geotiff driver supports 16bit, 32bit and 64bit floating point as well:

http://www.gdal.org/frmt_gtiff.html

A common usage are DEM files, which contain heigths as 16bit integer single band values.

This might help you on the topic http://www.iro.umontreal.ca/~mignotte/IFT6150/ComplementCours/UnderstandingRasters.pdf

And how to do it in python:

Python GDAL: Save array as raster with projection from other file

0

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.

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.