I wrote a script to rasterize different vector files the one after the other, during this process you have to input the cost per pixel for that specific layer. Eventually I am trying to achieve one raster file containing all the initial vector files which are now all represented by a different cost in a single raster file.
Now the problem is that when my inputted cost exceeds 255, that the cost in the raster file will just be 255. This is probably because of gdal using the Byte data-type. I tried to change it to GDAL_UInt16 to gat values up to 65535 but this is not working... When I output a file with the code below, I'll still get a file with values that range between 0 and 255.
Can anybody help me?
# -*- coding: utf-8 -*-
from osgeo import gdal, ogr, osr
import os
RASTERIZE_COST_FIELD = "__pixelkost__"
xmi = 0
xma = 0
ymi = 0
yma = 0
def rasterize(pixel_size=1):
# Ask for systempath towards main file
print ""
print "**************************************************"
print "** **"
print "** Please enter the systempath towards the **"
print "** shapefile which ranges over the total sur- **"
print "** face area of the used region. This will **"
print "** probably be the file with administrative **"
print "** lots. **"
print "** **"
print "**************************************************"
print ""
path = one()
# Ask for cost
cost1 = three()
# produce the outputpath
(shapefilefilepath, shapefilename) = os.path.split(path)
(shapefileshortname, extension) = os.path.splitext(shapefilename)
ofp = shapefilefilepath + "/" + shapefileshortname + ".img"
# Open the data source
orig_data_source = ogr.Open(path)
# Make a copy of the layer's data source because we'll need to
# modify its attributes table
source_ds = ogr.GetDriverByName("Memory").CopyDataSource(
orig_data_source, "")
source_layer = source_ds.GetLayer(0)
source_srs = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
xmi = x_min
ymi = y_min
xma = x_max
yma = y_max
# Create a field in the source layer to hold the features pixelkost
field_def = ogr.FieldDefn(RASTERIZE_COST_FIELD, ogr.OFTReal)
source_layer.CreateField(field_def)
source_layer_def = source_layer.GetLayerDefn()
field_index = source_layer_def.GetFieldIndex(RASTERIZE_COST_FIELD)
# Generate random values for the 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, cost1)
source_layer.SetFeature(feature)
# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(ofp, x_res,
y_res, 1, gdal.GDT_UInt16)
target_ds.SetGeoTransform((
x_min, pixel_size, 0,
y_max, 0, -pixel_size,
))
if source_srs:
# Make the target raster have the same projection as the source
target_ds.SetProjection(source_srs.ExportToWkt())
else:
# Source has no projection (needs GDAL >= 1.7.0 to work)
target_ds.SetProjection(ofp)
# Rasterize
err = gdal.RasterizeLayer(target_ds, [1], source_layer,
burn_values=[2509],
options=["ATTRIBUTE=%s" % RASTERIZE_COST_FIELD])
target_ds = None
Array1 = raster2array(ofp)
a = 1
print ""
print "**************************************************"
print "** **"
print "** Please enter, one by one, the systempath **"
print "** towards the shapefile you want to add to **"
print "** the set of rasterfiles. These will be **"
print "** merged to one final cost-file. **"
print "** **"
print "** If you want to stop adding shapefiles, type **"
print "** STOP **"
print "** **"
print "**************************************************"
print ""
while a != "STOP":
# Ask for systempath towards next file
path2 = one()
# Ask for cost
cost2 = three()
print ""
print "This process may take a while, please be patient."
# produce the outputpath
(shapefilefilepath, shapefilename) = os.path.split(path2)
(shapefileshortname, extension) = os.path.splitext(shapefilename)
ofp2 = shapefilefilepath + "/" + shapefileshortname + ".img"
# Open the data source
orig_data_source = ogr.Open(path2)
# Make a copy of the layer's data source because we'll need to
# modify its attributes table
source_ds = ogr.GetDriverByName("Memory").CopyDataSource(
orig_data_source, "")
source_layer = source_ds.GetLayer(0)
source_srs = source_layer.GetSpatialRef()
# Create a field in the source layer to hold the features pixelkost
field_def = ogr.FieldDefn(RASTERIZE_COST_FIELD, ogr.OFTReal)
source_layer.CreateField(field_def)
source_layer_def = source_layer.GetLayerDefn()
field_index = source_layer_def.GetFieldIndex(RASTERIZE_COST_FIELD)
for feature in source_layer:
feature.SetField(field_index, cost2)
source_layer.SetFeature(feature)
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(ofp2, x_res,
y_res, 1, gdal.GDT_UInt16)
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(ofp2)
err = gdal.RasterizeLayer(target_ds, [1], source_layer,
burn_values=[2509],
options=["ATTRIBUTE=%s" % RASTERIZE_COST_FIELD])
target_ds = None
Array2 = raster2array(ofp2)
length = len(Array1)
length2 = len(Array1[1])
l = 0
while l < length:
l2 = 0
while l2 < length2:
Array1[l,l2] += Array2[l,l2]
l2 += 1
l2 = 0
l += 1
while True:
try:
a = str(input("Press enter if you want to continue inputting shapefiles. Otherwise type STOP: "))
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 ""
return (Array1, path, ofp)
def one():
while True:
try:
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 ""
return syspath
def three():
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 ""
return cost
def raster2array(rasterfn):
raster = gdal.Open(rasterfn)
band = raster.GetRasterBand(1)
array = band.ReadAsArray()
return array
def array2raster(newRasterfn,rasterfn,array):
raster = gdal.Open(rasterfn)
geotransform = raster.GetGeoTransform()
originX = geotransform[0] # East/West location of Upper Left corner
originY = geotransform[3] # North/South location of Upper Left corner
pixelWidth = geotransform[1] # X pixel size
pixelHeight = geotransform[5] # Y pixel size
cols = array.shape[1]
rows = array.shape[0]
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, gdal.GDT_UInt16)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
if __name__ == '__main__':
A1, pa, imgpa = rasterize()
(shapefilefilepath, shapefilename) = os.path.split(pa)
ofp3 = shapefilefilepath + "/merge.img"
array2raster(ofp3,imgpa,A1)
print ""
print "**************************************************"
print ""
print "The merged rasterfile has been saved at:"
print ""
print str(ofp3)
print ""
print "**************************************************"
