python tips: array computation with numpy II

As previous post shows, for-loop calculation in array computations needs more processing time than vectorized calculation. There are some techniques for fasten calculation over arrays. Sample code shown below is one of them.

>>> from numpy import *
>>> 
>>> A = array([[1,2,3], [4,5,6], [7,8,9]])
>>> 
>>> def somefunc(x):
...     x1 = x
...     x2 = x - 3
...     return where(x < 7, x1, x2)
... 
>>> somefunc(A)
array([[1, 2, 3],
       [4, 5, 6],
       [4, 5, 6]])
>>> 

python tips: array computation with numpy

I compared performance of array computations. Listed below are results and codes of some experiments which introduced in the book named 'Python Scripting for Computational Science (3rd edition)'.

>>> from numpy import *
>>> import time
>>> 
>>> a = linspace(0, 1, 1E+07)
>>> t0 = time.clock()
>>> b = 3*a - 1
>>> 
>>> t1 = time.clock()
>>> 
>>> for i in xrange(a.size):
...     b[i] = 3*a[i] -1
... 
>>> t2 = time.clock()
>>> 
>>> print '3*a-1: %g sec, loop: %g sec' %(t1-t0, t2-t1)
3*a-1: 0.05 sec, loop: 15.46 sec

python tips: create shapefile with python

Long time has been passed since last post.

Here, I want to show the way to create Shapefile, which is standard data format for geometry data, using python.

Among several python tools such as geoscript, pyshp and pyproj which support operations of Shapefiles, I introduce how to use pyshp.

Pyshp is open source, provides shapefile library written all in pure python and distributed under MIT lisence.

You can download it as shapfile.py from  the official web site or just install via easy_install command:

# easy_install pyshp

However most of available operations of pyshp are introduced in the web site, I write some code using pyshp here just for my memo.

If you want to know more, I recommend looking over the web site and make some sample code according to tutorials.

#!/usr/bin/env python

import shapefile

def create_square(x=0., y=0., i=1.):
    return [[[x, y], [x+i, y], [x+i, y+i], [x, y+i], [x, y]]]

def define_geometry(w):
    # 1st record's geometry
    w.poly(parts=create_square(0, 0, 1))
    # 2nd record's geometry
    w.poly(parts=create_square(1, 1, 1))
    # 3rd record's geometry
    w.poly(parts=create_square(0.5, 0.5, 1))
    return w

def define_fields(w):
    w.field('FLD_1st')
    w.field('FLD_2nd', 'N', '10')
    return w

def fill_attributes(w):
    w.record('hoge', 'hogehogeho')
    w.record('mage', 'mage1')
    w.record('moge', 'mogemoge')
    return w

def create_proj_file(filename):
    # this function creates .proj file which defines projection.
    prj = open("%s.prj" % filename, "w")
    epsg = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
    prj.write(epsg)
    prj.close()

def main():
    filename = 'test'
    w = shapefile.Writer(shapefile.POLYGON)
    w = define_geometry(w)
    w = define_fields(w)
    w = fill_attributes(w)
    w.save('./' + filename)
    create_proj_file(filename)

if __name__ == "__main__":
    main()