How I implemented 2d affine transform in python

Although there is nothing special about it, but, strangely, I could not find it in GDAL-python (though it is very good library for working with raster and shape files, and I like it), also I looked at another respected library - shapely, but unfortunately, there is no capability to move, rotate or translate the shapes (but still there are some very handy functions like union, intersection, distance, area, ..). Another thing is that this can be easily done in Java using JTS, like this:
first I construct a polygon at the origin of the coordinate system and then determine the transformation needed to obtain this unit rectangle from the initial polygon (It can be defined uniquely using height of the polygon and the coordinates of the centers of the parallel edges). Then using the inverse transform we obtain the geometry object corresponding to the specification above (height, coordinates of the centers of the parallel sides). I drew these polygons using their wkt representation in OpenJump.

	import com.vividsolutions.jts.geom.Coordinate;
	import com.vividsolutions.jts.geom.Geometry;
	import com.vividsolutions.jts.geom.GeometryFactory;
	import com.vividsolutions.jts.geom.LinearRing;
	import com.vividsolutions.jts.geom.util.AffineTransformationBuilder;
	import com.vividsolutions.jts.geom.util.NoninvertibleTransformationException;
	import org.geotools.geometry.jts.JTSFactoryFinder;
	/**
	*
	* @author huziy
	*/
	public class TestCreateGeometry {
	public static void main(String[] args)
	throws NoninvertibleTransformationException {
	GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
	double width = 10.0;
	Coordinate c1 = new Coordinate(6,7);
	Coordinate c2 = new Coordinate(10,10);
	double height = c1.distance(c2);
	//vector perpendicular to c1c2
	Coordinate c = new Coordinate((c1.y - c2.y) / height * width / 2.0 , (c2.x - c1.x) / height * width / 2.0);
	Coordinate middle = new Coordinate(0.5 * (c1.x + c2.x), 0.5 * (c1.y + c2.y));
	Coordinate perpendicular = new Coordinate(middle.x + c.x, middle.y + c.y);
	Coordinate perpendDest = new Coordinate(0.5, 1.0); //impose it to be like that in the transformed plane
	//transform the points to the coordinates with the lower left at the corner of the crs
	//and edges parellel to the axes
	Coordinate c1Dest = new Coordinate(0, 0.5);
	Coordinate c2Dest = new Coordinate(1, 0.5);
	AffineTransformationBuilder builder = new AffineTransformationBuilder(c1, c2, perpendicular, c1Dest,
	c2Dest, perpendDest);
	//work from point1 and create linear ring of points
	Coordinate[] coords = new Coordinate[]{
	new Coordinate(0,0), new Coordinate(0,1), new Coordinate(1,1), new Coordinate(1,0), null
	};
	coords[4] = coords[0];
	LinearRing ring = geometryFactory.createLinearRing(coords);
	Geometry rect = geometryFactory.createPolygon(ring, null);
	rect = builder.getTransformation().getInverse().transform(rect);
	System.out.println("Perimeters");
	System.out.println(rect.getLength());
	System.out.println(2 * (height + width));
	System.out.println("Areas");
	System.out.println(rect.getArea());
	System.out.println(width * height);
	}
	}

view raw java_affine.java hosted with ❤ by GitHub

And still I decided to use python, since there exists an easy way of displaying results in matplotlib. Some people use descartes library for this purpose, but I don't feel like installing another dependency, since the thing with python is that if you change the OS, you have to reinstall them all. For me it is: matplotlib, numpy, scipy, gdal, cdat, shapely, netcdf4-python.... Well, I am getting distracted from the main point. So briefly, I feed my geometries to matplotlib as follows:

	polygon = loads(geom.ExportToWkt())
	boundary = polygon.exterior
	coords = np.zeros(( len(boundary.coords), 2))
	for i, the_coord in enumerate(boundary.coords):
	coords[i, 0], coords[i, 1] = basemap( the_coord[0], the_coord[1] )
	result.append(Polygon(coords, facecolor = 'none', edgecolor = edge_color, linewidth = linewidth))

view raw polygon_to_mpl.py hosted with ❤ by GitHub

So the transformation I did, can be initialized using 3 source and 3 corresponding destination points, that shoud not lie on the same line. Tested it, kind of works. Probably will optimize and improve it later, if necessary, but for now it looks like this:

	__author__="huziy"
	__date__ ="$Aug 10, 2011 3:49:16 PM$"
	from shapely.geometry import Point, Polygon
	import numpy as np
	class AffineTransform():
	def __init__(self):
	self.transMatrix = None
	self.shift = None
	pass
	def define_from_points(self, sourceToDest = {}):
	if len(sourceToDest) != 3:
	raise ValueError('only 2d transform is implemented')
	rhs = np.matrix(6 * [0.0]).transpose()
	mat = np.matrix(np.zeros((6,6)))
	sourcePoints = sourceToDest.keys()
	for i in xrange(0, 6, 2):
	source = sourcePoints[i / 2]
	dest = sourceToDest[source]
	# @type dest Point
	rhs[i] = dest.x
	rhs[i + 1] = dest.y
	for j in xrange(2):
	mat[i + j, 0 + 3 * j] = source.x
	mat[i + j, 1 + 3 * j] = source.y
	mat[i + j, 2 + 3 * j] = 1.0
	v = mat.getI() * rhs
	self.transMatrix = np.matrix(np.zeros((2,2)))
	self.transMatrix[0,0] = v[0,0]
	self.transMatrix[0,1] = v[1,0]
	self.transMatrix[1,0] = v[3,0]
	self.transMatrix[1,1] = v[4,0]
	self.shift = np.matrix([v[2,0], v[5,0]]).transpose()
	def transform(self, geometry):
	if type(geometry) == type(Point()):
	v = np.matrix([geometry.x, geometry.y]).transpose()
	res = self.transMatrix * v + self.shift
	return Point(res[0,0], res[1,0])
	if type(geometry) == type(tuple([])):
	return self.transform(Point(geometry))
	elif type(geometry) == type(Polygon()):
	# @type geometry Polygon
	ps = geometry.exterior.coords
	newPs = []
	for p in ps:
	thePoint = self.transform(p)
	newPs.extend(thePoint.coords)
	return Polygon(newPs)
	else:
	raise ValueError('Unknown geometry type: {0}'.format(geometry) )
	def test_translation():
	aft = AffineTransform()
	points = {}
	p1 = Point(0,0)
	p2 = Point(1,1)
	p3 = Point(0,1)
	points[p1] = Point(1,1)
	points[p2] = Point(2,2)
	points[p3] = Point(1,2)
	aft.define_from_points(points)
	p1 = aft.transform(Point(10,10))
	assert Point(11.0,11.0).wkb == p1.wkb
	def test_rotation():
	aft = AffineTransform()
	points = {}
	p1 = Point(0,0)
	p2 = Point(1,1)
	p3 = Point(0,1)
	points[p1] = Point(0,0)
	points[p2] = Point(-1,-1)
	points[p3] = Point(0,-1)
	aft.define_from_points(points)
	assert aft.transform(Point(10,10)).wkb == Point(-10, -10).wkb
	def test_polygon():
	aft = AffineTransform()
	points = {}
	p1 = Point(0,0)
	p2 = Point(1,1)
	p3 = Point(0,1)
	points[p1] = Point(0 + 10,0 + 15)
	points[p2] = Point(-1 + 10,-1 + 15)
	points[p3] = Point(0 + 10,-1 + 15)
	pts = []
	pts.extend(p1.coords)
	pts.extend(p2.coords)
	pts.extend(p3.coords)
	pts.extend(p1.coords)
	polygon = Polygon([(-1, -1), (1,-1), (1,1), (-1,-1)])
	aft.define_from_points(points)
	print polygon.wkt
	p1 = aft.transform(polygon)
	print p1.wkt
	def test():
	test_translation()
	test_rotation()
	test_polygon()
	if __name__ == "__main__":
	test()
	print "Hello World"

view raw shapely_affine.py hosted with ❤ by GitHub

Basically, I just solve a system of linear equations with respect to parameters of the transformation, and then use the parameters to transform other points.