So far I have developed the below code but for some reason it only works for say one class of land-use instead of two. Would you please tell me where I have gone wrong?

import os, numpy, osgeo.gdal, math, random, sys
from osgeo.gdal import *
import  osgeo.gdalconst, osgeo.ogr
driver = osgeo.gdal.GetDriverByName('RST')
ds = osgeo.gdal.Open('landii.rst')
cols = ds.RasterXSize
rows = ds.RasterYSize
print 'no of cols is:', cols
print 'no of rows is:', rows
print 'Raster has %s bands' % ds.RasterCount
geotransform = ds.GetGeoTransform()

out = driver.Create ('chandal.rst', cols, rows, 1, osgeo.gdal.GDT_Byte)
driver = ds.GetGeoTransform
raster = numpy.zeros((rows, cols), dtype=numpy.int32)

for i in range(rows):
    for j in range(cols):
        coll = 0
        raster[i][j] = coll
perc = [0.65, 0.45] # perc is the percentage of the pixels to be allocated to land use 1 and two respectively

nclasses = 2
for ii in range(nclasses):
    l = perc[ii]
number = int(cols * rows * l)
while (number > 0):
    rand1 = random.randrange(rows)
    rand2 = random.randrange(cols)
    k = raster[rand1][rand2]
    if k == 0:
        raster[rand1][rand2] = ii 
        number = number - 1


out = None
  • you have all the spaces. Is this just because something went wrong copying it in here or is this actually your code. Maybe you could updated it, before we start working on the code. – ustroetz Apr 16 '15 at 10:42
  • Sorry, I forgot to mention that the spaces here are mistakes, – Rassoul Apr 16 '15 at 12:32
  • Maybe you could update the code without the spaces, this way it would be easier to work with it. – ustroetz Apr 16 '15 at 12:34
  • please forget about the lines spaces – Rassoul Apr 16 '15 at 12:44

Your best bet to get a weighted choice is to use numpy.random.choice which allows you to specify the sample set, and a weight for each sample. Note the probability must sum to exactly 1.

raster = numpy.random.choice([0, 1], size=(rows, cols), p=[0.65, 0.35])

Also, a quick note: your comment says the probability that the landuse should be allocated to 1 or 2, but your data type is a binary raster which can only ever have the values of 0 or 1. Instead you should use an integer type if you want to keep the values 1 and 2; gdal.GDT_Int32

Also, if you're looking at expanding on your code above a couple of code issues:

  • There is no need to pre-set the zero raster values to zero - it's already done when the array is initialised with numpy.zeros.
  • With a consistent input raster the value of number will always be the same - effectively cols * rows * 0.45. The loop over the perc finishes before you assign the value to number.
  • You can index a 2D numpy array with raster[row, col]
  • Finally rather than using GDAL directly I'd take a look at the rasterio library which is a much more pythonic wrapper around GDAL.
| improve this answer | |

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.