## Three good StackOverflow answers on Python topics

Here are three thorough answers on Python topics from StackOverflow.

## 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```

## 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() srs.ImportFromEPSG(25832) 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.

## Using Cython, psyco and other tricks to make Python run faster

This is a good post showing how “perrygeo” uses Cython to speed up a pure Python function. The comments also discuss alternatives and options:

• Using psyco instead of Cython
• Using optimization flags for gcc
• Using xxx instead of math.xxx

http://www.perrygeo.net/wordpress/?p=116

## 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```

## Transforming geographical and projected coordinates in Python

pyproj is a Python interface to the PROJ.4 (http://trac.osgeo.org/proj/) functions. It allows you to transform coordinates between coordinate systems.

## Install pyproj using pip

Assuming you’re using pip to install Python packages:

`\$ pip install pyproj`

## Python module for processing a file line-by-line

Note: Since writing this post, I’ve learned about the fileinput module, which turns most of the following into a oneliner:

```import fileinput for line in fileinput.input(): process(line)```

## 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 gis.stackexchange.com

gis.stackexchange.com 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
http://gis.stackexchange.com/questions/19401/projecting-shapefile-with-transformation-using-ogr-in-python

Will add more as I find them…