Welcome, guest | Sign In | My Account | Store | Cart

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
  1
  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.

Created by Gerry Wiener on Thu, 12 May 2005 (PSF)
Python recipes (4591)
Gerry Wiener's recipes (4)

Required Modules

Other Information and Tasks