Creating visually random images

How complicated does a mathematical function, pseudorandom(x), have to be to create something that seems random, and how do you check whether it seems random? Looking at a long list of numbers is not a good way, because our ability to look at long lists of numbers is very limited. Our ability to look at images is much better.

So, to inspect whether a list of numbers is seemingly random, a fun way is to create an image using the numbers for each pixel, and simply look at it.

Given a number x in the series {0, 1, 2, …, huge n}, let’s create an image that plots a random function that’s implemented using the sin() function and a 3rd degree polynomial function:

color_range = 2**8
a = -3.0
b = -5.0
c = 13.0
def pseudo_randombit(x):
	color255 = math.sin(a*x**3 + b*x**2 + c*x) % color_range
	# make black/white
	bit = color255 / 127
	return bit

Continue reading “Creating visually random images”

Incrementing part of an array in numpy

given a (2D) array a, and a smaller (2D) array b, how can you add the smaller array to the bigger array at some offset in numpy?

>>> a = zeros((4,4), dtype=int)
>>> a
array([[0, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 0, 0]])
>>> b = ones((2,2), dtype=int)
>>> b
array([[1, 1],
       [1, 1]])
>>> a[1:3,1:3] += b
>>> a
array([[0, 0, 0, 0],
       [0, 1, 1, 0],
       [0, 1, 1, 0],
       [0, 0, 0, 0]])

Hello world of raster creation with GDAL and Python

Mostly so I myself can remember how to do it, here is how to create a random geotiff with GDAL in Python

Note: the width and height are given in opposite order in the GDAL raster and numpy arrays!

import osr
import numpy
import gdal
import math
width = 4000
height = 3000
format = "GTiff"
driver = gdal.GetDriverByName( format )
dst_ds = driver.Create( "test.tiff", width, height, 1, gdal.GDT_Byte )
dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )
srs = osr.SpatialReference()
dst_ds.SetProjection( srs.ExportToWkt() )
raster = numpy.zeros( (height, width), dtype=numpy.uint32 )
color_range = 2**8
seed = math.pi**10
for i in range(height):
	for j in range(width):
		color = (seed*i*j) % color_range
		raster[i][j] = color
dst_ds.GetRasterBand(1).WriteArray( raster )

It’s kind of slow, so perhaps the operation can be speeded up somehow? The result looks kind of nice though (image created with width and height both 4000):

So, not completely random.

Quick benchmarking of Python code with timeit

How long does [insert code snippet] take in Python? Find out with the timeit module, in Python since 2.3.

Below I time how long it takes to make a list with a 1000 zeros, [0,0,0,…], using different Python code snippets. Notice that you can time snippets directly on the commandline, without writing a script file (.py) or running the Python interpreter interactively.

Benchmarks ordered by running time

Using numpy:

$ python -m timeit 'import numpy' 'a = numpy.zeros(shape=(1,1000))'
100000 loops, best of 3: 2.38 usec per loop

Continue reading “Quick benchmarking of Python code with timeit”

GDAL and OGR Python tips

Online documentation for GDAL/OGR Python is sparse. Here I show some recommended ways of learning more about GDAL/OGR in Python.

Using Python interpreter

You can learn about GDAL and OGR from inside the Python interpreter.

Start python interpreter:

$ python

Import the modules:

from osgeo import gdal,ogr,osr

Learn about the modules using ‘help’ and ‘dir’ in Python. These built-in functions work for any type of object (module, class, functions etc):

help(ogr) # display help for the ogr module
dir(ogr) # what's contained in ogr module?
help(ogr.Geometry) # display help for ogr.Geometry class
dir(ogr.Geometry) # show contents of an object, like functions on a class

Using is a Q&A site about GIS. Many people here know about GDAL and ORG, also about the Python bindings.

Some good questions I’ve found:

Transformation using OGR in Python

Will add more as I find them…