AGU Fall Meeting 2011 and around that time

View from the airport (SFO)

Line to buy BART tickets

My hotel room

Moscone West hall

Poster Hall

My poster

Alkatraz

Golden Gate bridge

find

Linux: search for a folder with name afb:

find ./ -name afb -type d

view raw find_ex1.sh hosted with ❤ by GitHub

To which cell the point belongs

Recently I was relating measurement points and grid cells, using shapely's Polygon.intersect method. And was lucky enough to bump into the situation where the point, which corresponds to a real observation station, was not inside any grid cell. Though it was inside the simulation domain judging by its coordinates. Fortunately we got Polygon.buffer method, which allowed me to extend the polygons, in order to make sure that every point inside the domain will fall into some cell. The length of the buffer was selected as 0.01 * (cell edge).

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.

Adding text files

If you have many files with text information, it is easier to look for what you need, if
all the info is assembled into a single file. To achieve this in Bash, you can do the following:

cat $(ls) > dictionary.txt

view raw merge_txt_files.sh hosted with ❤ by GitHub

this will put all files from the current directory into one file.
More complicated manipulations, such as sorting and elimination of duplicates,
described here.

Bash: kill processes with a given name pattern

As we know the kill command in linux permits to kill a process using its id. kill -9 tries to kill your process even harder. But what if I have many similar processes? We can use the cut function that I recently saw here as shown below:

	for id in $(ps ax | grep bootstrap.py | cut -d" " -f1)
	do
	kill -9 $id
	done

view raw kill.sh hosted with ❤ by GitHub

What it does is: searches the list of processes for lines containing bootstrap.py, then extracts the process id from each line and kills the process.

C++: Allocating 2D array using new

This link explains how to allocate 2d arrays with variable dimensions:
the link. I find that it is a very useful and clear explanation.

NCO

NCO is a very powerfull tool for editing netcdf files. In this post I will be writing some useful commands.

1. Delete a dimension of size 1 (the dimension name is level):

ncwa -a level directions_qc_amno.nc directions_qc_amno1.nc

2. Copy data from one to the other netcdf file, if possible recognize common dimensions:

ncks file1.nc file2.nc

3. Rename a dimension

ncrename -d rlat,latitude directions_qc_amno.nc

4. permute dimensions(make longitude first and the latitude second):

ncpdq -a longitude,latitude infocell.nc infocell.nc 

5. Rename the variable lat to latitude:

ncrename -v lat,latitude infocell.nc

6. Take difference of variable x in 2 files:

ncbo -v x --op_typ='-' infocell.nc infocell_1.nc infocell_diff.nc

7. Operation on variables in a file. Here the values of the variable are changed and the attributes added.

ncap -O -s "air=(air-512.81)/0.01" filename.nc temp.nc
ncap -O -s "air=short(air);air@add_offset=512.81;air@scale_factor=0.01" temp.nc filename.nc

taken from http://jisao.washington.edu/data/nco/
8. Concantenating netcdf files along the record dimension.

ncrcat runoff_file_name_pattern total_runoff_af.nc

9. Subsetting by index and by coordinate values:

ncks -d lon,0,4 -d lat,2,9,2 in.nc out.nc
ncks -d lon,10.,20. -d lat,80.,90. in.nc out.nc

Used data from:

• http://linux.die.net/man/1/ncbo
• http://nco.sourceforge.net/#RTFM