This recipe illustrates how to extract a low resolution grid from a high resolution grid using the Numeric package. It presents an inefficient but straightforward solution and then a more efficient solution that employs the functionality in Numeric.
Python, 129 lines
Download
Copy to clipboard1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 | #!/usr/bin/env python
"""Replace a grid with one that is more coarse"""
import Numeric
def grid_coarse(in_grid, in_ny, in_nx, out_ny, out_nx, factor):
"""Replace a grid with one that is more coarse
Inputs:
in_ny -- number of input grid rows
in_nx -- number of input grid columns
out_ny -- number of output grid rows
out_nx -- number of output grid columns
factor -- input grid box length divided by output grid box length
"""
out_grid = Numeric.zeros((out_ny, out_nx), in_grid.typecode())
for j in xrange(0, out_ny):
xjfine = float(j) * factor
jfine = int (float(j) * factor)
dj = xjfine - float(jfine)
jp1 = min(in_ny-1, jfine+1)
jm1 = max(0, jfine-1)
for i in xrange(0, out_nx):
xifine = float(i) * factor
ifine = int (float(i) * factor)
di = xifine - float(ifine)
ip1 = min(in_nx-1, ifine+1)
xval = (1.-di) * (1.-dj) * in_grid[jfine, ifine] + (1.-di) * dj * in_grid[jp1, ifine] + di * (1.-dj) * in_grid[jfine, ip1] + di * dj * in_grid[jp1, ip1]
out_grid[j, i] = xval
return out_grid
class GridCoarse(object):
"""Replace a grid with one that is more coarse"""
def __init__(self, in_ny, in_nx, out_ny, out_nx, factor):
"""Constructor
Inputs:
in_ny -- number of input grid rows
in_nx -- number of input grid columns
out_ny -- number of output grid rows
out_nx -- number of output grid columns
factor -- input grid box length divided by output grid box length
"""
self.in_ny = in_ny
self.in_nx = in_nx
self.out_ny = out_ny
self.out_nx = out_nx
self.factor = factor
self.jfine = Numeric.zeros(out_ny)
self.dj = Numeric.zeros(out_ny, Numeric.Float64)
self.jp1 = Numeric.zeros(out_ny)
self.ifine = Numeric.zeros(out_nx)
self.di = Numeric.zeros(out_nx, Numeric.Float64)
self.ip1 = Numeric.zeros(out_nx)
for j in xrange(0, out_ny):
xjfine = float(j) * self.factor
self.jfine[j] = int(float(j) * self.factor)
self.dj[j] = xjfine - float(self.jfine[j])
self.jp1[j] = min(self.in_ny-1, self.jfine[j]+1)
for i in xrange(0, out_nx):
xifine = float(i) * self.factor
self.ifine[i] = int(float(i) * self.factor)
self.di[i] = xifine - float(self.ifine[i])
self.ip1[i] = min (self.in_nx-1, self.ifine[i]+1)
self.d1 = Numeric.outerproduct(1.0 - self.dj, 1.0 - self.di)
self.d2 = Numeric.outerproduct(self.dj, 1.0 - self.di)
self.d3 = Numeric.outerproduct(1.0 - self.dj, self.di)
self.d4 = Numeric.outerproduct(self.dj, self.di)
def convert(self, in_grid):
g1 = Numeric.take(Numeric.take(in_grid, self.jfine), self.ifine, 1)
g2 = Numeric.take(Numeric.take(in_grid, self.jp1), self.ifine, 1)
g3 = Numeric.take(Numeric.take(in_grid, self.jfine), self.ip1, 1)
g4 = Numeric.take(Numeric.take(in_grid, self.jp1), self.ip1, 1)
out_grid = self.d1 * g1 + self.d2 * g2 + self.d3 * g3 + self.d4 * g4
return out_grid.astype(in_grid.typecode())
def test1():
in_ny = 500
in_nx = 500
out_ny = 250
out_nx = 250
factor = 2.0
a = Numeric.arange(500*500)
in_grid = a.resize((500, 500))
gc = grid_coarse(in_grid, in_ny, in_nx, out_ny, out_nx, factor)
def test2():
in_ny = 500
in_nx = 500
out_ny = 250
out_nx = 250
factor = 2.0
a = Numeric.arange(500*500)
in_grid = a.resize((500, 500))
gc = GridCoarse(in_ny, in_nx, out_ny, out_nx, factor)
out_grid = gc.convert(in_grid)
if __name__ == "__main__":
import timeit
t = timeit.Timer("test1()", "from __main__ import test1")
print "test1(): ",
print t.timeit(number=10)
t = timeit.Timer("test2()", "from __main__ import test2")
print "test2(): ",
print t.timeit(number=10)
|
This recipe illustrates how functionality in the Numeric package can be applied in cases where one would normally iterate through a rectangular subgrid. The initial solution is presented in grid_coarse(). This solution is quite efficient in Fortran or C but is not an effective solution when used with the Numeric package. The convert method in GridCoarse avoids the inner loop and is useful for converting a number of 2-d grids all having the same grid dimensions. This is just what is called for when converting the horizontal planes in a 3d volume.
The original Fortran version for grid_coarse() was provided by the Forecast System Laboratory at NOAA.
Tags: graphics
