[Image-SIG] Raster Calculation: Script Error

[Fredrik Lundh]

On Fri, Apr 3, 2009 at 4:59 AM, santosh panda <santoshh19 at yahoo.co.in> wrote:

> I am using Landsat single band image as input and I want to apply the
> equation shown below in the script and write the output to a new raster
> file.
>
> 1. The below pasted script works fine without 'cos' term but it fails with
> 'cos' term even if I imported python math module. Should I import any other
> module so that the script will recognize 'cos' term?
>
> 2. I want my output data to be in same projection and pixel size as my input
> data. What command do I have to use to maintain the projection details and
> pixel size?
>
> 3. My input data is 8 bit and I want my output data also in 8 bit. How can I
> maintain the data mode?
>
> script:
> from PIL import Image, ImageMath
> import os, math
> os.curdir = 'C:/temp'
> im = Image.open("testraster.tif")
> out = ImageMath.eval("((3.14*(((0.11*im) + lmin)*25))/(70*cos(50*
> 0.0174)))", im=im, lmax= 30, lmin= 0.37)
> out.save("outRaster.tif")

The ImageMath eval function only supports a small number of functions.
 You can either calculate that part of the expression outside the call
and pass it in as a parameter, or pass in cos() itself as a callable
parameter:

    out = ImageMath.eval(..., cos=math.cos)

The eval function uses I or F as intermediate modes; to convert back
to a smaller mode, you can either use the convert function as
described here:

    http://effbot.org/imagingbook/imagemath.htm#built-in-functions

or do the conversion before saving the image:

    out = ImageMath.eval(...)
    out = out.convert("L")
    out.save("...")

>From you example above, it looks like the expression is constant for
the entire image; in that case, using the "point()" method on the
Image object is probably more efficient (at least as long as you're
working on 8-bit images).  You can either build a lookup table, or
pass in a lambda or other callable; point will then calculate the
value once for each possible pixel value, and translate the image in
one pass over all pixels: e.g.

   def function(pixel):
      return math.pi * (((0.11*pixel) + lmin)*25))/(70*math.cos(50*> 0.0174)))

   out = im.point(function)

This will call the function once for each possible value, and then do
the mapping.  If you're processing a large number of images, you can
even precompute the table:

   lut = map(function, range(256))

and then do "out = im.point(lut)" for each image.

</F>