FFT using pylab

from pylab import *

data,timestep = 10*sin(arange(1000)*5/2.0/pi)+5*sin(arange(1000)*10/2.0/pi),0.01
han = hanning(len(data))
handata = han*data
n = len(data)*50
FFT = abs(fftshift(fft(handata,n)))
freq = fftshift(fftfreq(n,timestep))
plot(freq,2*FFT/sum(han))

RoiView - Explore InSAR data and more

RoiView was originaly designed to replace DGX/MDX software which require bureaucracy to obtain. This software was used in order to view ROI-PAC interferometry processing results but producing an image from file was not so strait forward. In time support was added for ERDAS ER Mapper header files and for user provided parameters with or without header files.

Current RoiView version (0.75) is available for download at
See also RoiView's help files

Screenshot:

Using hillshade image as intensity (improved matplotlib shade)

Matplotlib module enables hillshade method (shade) using a LightSource class (v 0.99).
The problem is it uses the data itself as intensity and data. It is very useful for viewing a DEM but sometimes you would like the DEM as intensity underlying some other data. Another problem is that the shade method is producing a very light colored image sometimes even white where intensity is high.
I used as an example a DEM derived from SRTM v4 data acquired at the International  Centre for Tropical  Agriculture (CIAT - http://srtm.csi.cgiar.org) the hillshade production was made using LightSource class with azimuth of 165 deg. and altitude of 45 deg.)

DEM - gist_earth color scheme

Hill-shade (azdeg-165,altdeg-45)

matplotlib shade method

My shade method

The difference in the shading colors derived from the method used to produce it. While the matplotlib method uses "hard light" method I use a "soft light" method. the matplotlib is converting the RGB colors to HSV and then calculate the new saturation and value according to the intensity. I use a formula based on the description of ImageMagick's pegtop_light.which is much faster as it is a single formula. Another advantage is the option to use a separate layer as the intensity and another as the data used for colors.
The modified functions are hillshade and set_shade as follows:

#!/bin/env python
from pylab import *
def set_shade(a,intensity=None,cmap=cm.jet,scale=10.0,azdeg=165.0,altdeg=45.0):
''' sets shading for data array based on intensity layer
  or the data's value itself.
inputs:
  a - a 2-d array or masked array
  intensity - a 2-d array of same size as a (no chack on that)
                    representing the intensity layer. if none is given
                    the data itself is used after getting the hillshade values
                    see hillshade for more details.
  cmap - a colormap (e.g matplotlib.colors.LinearSegmentedColormap
              instance)
  scale,azdeg,altdeg - parameters for hilshade function see there for
              more details
output:
  rgb - an rgb set of the Pegtop soft light composition of the data and 
           intensity can be used as input for imshow()
based on ImageMagick's Pegtop_light:
http://www.imagemagick.org/Usage/compose/#pegtoplight'''
  if intensity is None:
# hilshading the data
    intensity = hillshade(a,scale=10.0,azdeg=165.0,altdeg=45.0)
  else:
# or normalize the intensity
    intensity = (intensity - intensity.min())/(intensity.max() - intensity.min())
# get rgb of normalized data based on cmap
  rgb = cmap((a-a.min())/float(a.max()-a.min()))[:,:,:3]
# form an rgb eqvivalent of intensity
  d = intensity.repeat(3).reshape(rgb.shape)
# simulate illumination based on pegtop algorithm.
  rgb = 2*d*rgb+(rgb**2)*(1-2*d)
  return rgb

def hillshade(data,scale=10.0,azdeg=165.0,altdeg=45.0):
  ''' convert data to hillshade based on matplotlib.colors.LightSource class.
    input:
         data - a 2-d array of data
         scale - scaling value of the data. higher number = lower gradient
         azdeg - where the light comes from: 0 south ; 90 east ; 180 north ;
                      270 west
         altdeg - where the light comes from: 0 horison ; 90 zenith
    output: a 2-d array of normalized hilshade
'''
  # convert alt, az to radians
  az = azdeg*pi/180.0
  alt = altdeg*pi/180.0
  # gradient in x and y directions
  dx, dy = gradient(data/float(scale))
  slope = 0.5*pi - arctan(hypot(dx, dy))
  aspect = arctan2(dx, dy)
  intensity = sin(alt)*sin(slope) + cos(alt)*cos(slope)*cos(-az - aspect - 0.5*pi)
  intensity = (intensity - intensity.min())/(intensity.max() - intensity.min())
  return intensity

Example of use:
One can save the code to a file named say: shading.py
now say we have a 4 byte float DEM data in a 560 lines 420 samples binary file. in a python code:

from pylab import *
from shading import set_shade
from shading import hillshade
dem = fromfile('DEM.dem',dtype=float32).reshape(560,420)
rgb = set_shade(dem,cmap=cm.gist_earth)
imshow(rgb)

will produce the "my shade method" image as above.

say we have a data to be plot using the DEM data as intensity:

replace the line before last with:
rgb = set_shade(data,intensity=hillshade(dem),cmap=cm.gist_earth)