Archive Page 2

install gdal on Mac OS X


Mac OS X is built on Unix BSD and thus offers terminal and bash for good old scripting.

To add gdal and its executables (gdal_translate, gdal_merge.py etc.), go first downloading the framework prepared by KyngChaos

http://www.kyngchaos.com/software:frameworks

and download Gdal Complete.dmg.

The adaptation done for Mac OS X  is just perfect! To make it run:

mount the dmg, then double click on the installation package. Once the installation is done, you must set the PATH variable for bash.

Launch Terminal, then type

export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH

This command adds /Library/Frameworks/GDAL.framework/Programs to the search PATH.

If you want to save this setting, edit the hidden file .bash_profile in your home directory, and add the above command line; then save. For beginners: the dot ‘.’ before bash_profile corresponds to a hidden file. To immediately set the path, type

source .bash_profile

Anyway, .bash_profile will be read again next time you run a terminal.

Enjoy!

Advertisements

interpolating colors

I’m starting a little project: making a python function to apply a color scheme to a gray level image.

For that we need: to read the image (job done by gdal), read a color scheme, get image min an max excluding no data values, interpolate the color scheme, and last but not least, transform the single band gray level into an rgb image.

For today, let’s start with the color interpolation function. Color can be interpolated in rgb, hsv, etc. Let’s write something for rgb, easily portable to other color space.

The idea is to consider separately the 3 color axis and interpolate along them:

def interpolateColor(rgbMin, rgbMax, maxDepth, depth):
    if depth <= 0:
        return rgbMin
    if depth >= maxDepth:
        return rgbMax

    rgb=rgbMin
    for iTriplet in range(3):
        minVal=rgbMin[iTriplet]
        maxVal=rgbMax[iTriplet]
        rgb[iTriplet] = int(minVal + (maxVal-minVal)*(depth/float(maxDepth)))

    return rgb

def test():
    rgbMin=[128, 0 ,34]
    rgbMax=[255, 56, 0]

    niter=10
    for ii in range(niter):
       print ii, interpolateColor(rgbMin, rgbMax, niter, ii)

To test it:

import interpolateColor
interpolateColor.test()

You should get the 9 intermediate colors.

python-gdal: simple file processing

Here is a simple example of image processing with python and gdal.

Say we want to make a mask from an image, 0 if data<10, 1 else.
The strategy is the following:

  • open the input file
  • create an output file, geotiff, byte
  • read the input file line by line
  • test values and write to the output

First, let’s import some libraries. I use the following snippet in any of my codes, so it must (may?) work with most of the gdal/numpy libraries settings:

try:
    from osgeo import gdal
    from osgeo.gdalconst import *
    gdal.TermProgress = gdal.TermProgress_nocb
except ImportError:
    import gdal
    from gdalconst import *
try:
    import numpy as N
    N.arrayrange = N.arange
except ImportError:
    import Numeric as N
try:
    from osgeo import gdal_array as gdalnumeric
except ImportError:
    import gdalnumeric

Now let’s define file names, and the threshold value. You can read something else than a geotiff file in input, gdal should do the job for you.

file='/data/demo/input.tif'
outfile='/data/demo/output.tif'
value=1

The lines below open the input file, and get from it its width and height.

fid = gdal.Open(file, GA_ReadOnly)
ns = fid.RasterXSize
nl = fid.RasterYSize

The output file is created, with the same dimensions, 1 single band, its type is set to Byte (GDT_Byte), format forced to geotiff, options are left empty. Note that options is an array of strings, we’ll see in another post what can be done with it. The output file get same projection and geotransformation as the input. Geo transformation is an array of float indicating the upper left corner coordinates and pixel size.

outDrv = gdal.GetDriverByName('GTiff')
options=[]
outDs  = outDrv.Create(outfile, ns, nl, 1, GDT_Byte, options)
outDs.SetProjection(fid.GetProjection())
outDs.SetGeoTransform(fid.GetGeoTransform())

Preparation is done, it is now enough to parse the input file line by line.
At each iteration the output array is set to 0 (default value), places where data is found greater than the given threshold being set to 1.

Note that ReadAsArray returns a 2-D array. N.ravel transforms it in a single array on which operations like dataOut[data>value]=1 can be applied. To write this array to a file you MUST give it 2 dimensions, which is achieved with dataOut.shape=(1,-1).

datazeros=N.zeros(ns)
for il in range(nl):
data = N.ravel(fid.GetRasterBand(1).ReadAsArray(0, il, ns, 1))
dataOut = datazeros
dataOut[data>value]=1
dataOut.shape=(1,-1)
outDs.GetRasterBand(1).WriteArray(N.array(dataOut), 0, il)

Job done!

Time machine freezed with a fire wire cable

Although I prefer my DYI backup system, I also tried Time Machine on my Mac Mini, on Leopard.

Problem: After making some save sessions ok, Time Machine was in ‘update’ status for hours without saving anything, forcing me to reboot the machine (could not even kill the app cleanly)!

I changed the disk format to GUID as suggest in http://support.apple.com/kb/TS1550?viewlocale=en_US, still had the problem.

I finally found it was due to the fire wire cable. I changed to USB2 connection: no problem anymore!

No clue if it was the electronic on board the external drive (no-name hardware) or compatibility or bug in Leopard, anyway changing the cable solved everything.

May worth to try.

Faster loop in Python

If found this today: I had a loop for scanning an image line by line:

for il in range(ystart, yend+1):

data = N.ravel(fid.GetRasterBand(band+1).ReadAsArray(xstart, il, ns,1))
wts = (data != 0)

if wts.any():

index = N.arange(ns*nl)
jpos  = index / ns
ipos  = index - jpos * ns
jmin = jpos.min()
jmax = jpos.max()
imin = ipos.min()
imax = ipos.max()
etc...

The loop was pretty slow.  I thought it was coming from computation on the arrays in the 'if' statement. Actually, by commenting those lines, I realized it was due to the creation of the index array:
index = N.arange(ns*nl)
To (significantly) speed up the loop, I simply had to create first a reference indexRef array to avoid creating it in the loop.

indexRef = N.arange(ns*nl)
for il in range(ystart, yend+1):

data = N.ravel(fid.GetRasterBand(band+1).ReadAsArray(xstart, il, ns,1))
wts = (data != 0)

if wts.any():

index = indexRef[wts]
jpos  = index / ns
etc...

Making an image layer stack for time series processing

One often need a layer stack image, at least for displaying time series of values for a pixel. Have you seen that Envi make_layerstack function provided with the GUI sorts the file name in reverse order? Quite irritating. If you’ve got gdal, simply make the layer stack in one command lines.

gdal_merge.py -of gtiff -separate -o layerstackname.tif file1 file2 file3 file4

I use to make layer stack of hundreds images for processing time series, which I got in two command lines:

list=($(ls *_ndvi.tif))

gdal_merge.py -of gtiff -separate -o timeSeries.tif ${list[@]}

The list variable is actually an array (the two external parenthesis). You can be more selective in the array definition by playing around with the ls command or adding some grep command.

Is a script running?

On a *nix system, it is often desirable to ensure that a script is running as a singleton (I borrow the C++ term). For example, if a process is scheduled to run every 10min and for some reason last for 11, for sure you don’t want it to be run again.

There are some tricks, like locking a socket and checking its status before running, but I guess it impose quite a lot of programmation, maybe a bit difficult for a shell script. A simpler way is to trigger the script from another one (let’s call it singleton) that is charged to check if the script (cunningly named myCommand) is already running or not, like:

singleton myCommand

What should do singleton? Well, there are two cases: either myCommand is an executable directly run by the system or it is a script. If it is an executable, then pidof let’s you check if it is running:

pidof -s -x myComnand

Now, if myCommand is a script, starting with a shabang (#!/bin/bash or #!/bin/env python, etc) pidof won’t find anything since the program actually running is bash or python and myCommand is only a parameter. Check it with ps -eaf | grep myCommand and you’ll see:

/bin/bash myPath/myCommand.sh

To check if myCommand is executed, the simplest is to list processes, grep lines with myCommand.sh. Don’t forget to add an inverse grep on grep itself (grep -v grep):

result=$(ps -eaf | grep -v grep | grep myCommand)

if [ -z “$result” ]; then

myCommand

fi

Of course, you can be less restrictive when selecting processes from ps. In this case, singleton won’t run anything if you are editing the script (example vi myCommand.sh) or simpling displaying it (e.g. cat myCommand.sh). It is up to you to set the restriction.