gdal_calc.py– Raster Calculator using Numby Functions

gdal_calc is a python script that makes it easy to do band math and logical operations with gdal/numby. This ranges from basic arithemtic operations to logical functions. And while gdal_calc.py has been around since 2008, it is but is a recent and revelatory discovery for me. I had just today resigned myself to properly learning python so as to use the gdal bindings. But, my laziness may continue unabated. Not learning Python. Now that is lazy.

Anyway, want to do band math? It’s as simple as follows (straight from the python script itself):

# example 1 - add two files together
gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"

# example 2 - average of two layers
gdal_calc.py -A input.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2"

# example 3 - set values of zero and below to null
gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0

So, I had little modification to average my set of images from PovRay:

gdal_calc.py -A input.tif -B input2.tif -C input3.tif -D input4.tif -E input5.tif --outfile=result.tif --calc="(A+B+C+D+E)/5"

I might write some bash wrappers to create built in functions for basic things I’ll do repeatedly, like average a set of images based on a directory listing. Of course, a little voice inside says “You should just do that in Python.” We’ll see.

GDAL, MrSid, and nearblack

Translating MrSid lossy compressed files into uncompressed imagery has its drawbacks, including licensing and artifacts.
Old versions of fwtools, which includes the GDAL utilities (and more), were compiled with a license that allowed for the translation of MrSid into e.g. Erdas Imagine images or GeoTiff. The licensing changed on that library, so FWTools and MS4W don’t do this anymore.

If you have a compiled version of the GDAL utilities with a legal license for converting, you still run into the issue of artifacts, e.g.

Nearly, but not quite black pixels at the edge of MrSid lossy compressed image

GDAL’s nearblack is a great utility for correcting this problem. If we wanted to batch this up on a windows system, it might look something like this:

FOR %f IN (*.sid) DO nearblack -o %~nf.img %f

Landscape Position Continued– absolutely relative position calculation <– Pics!

Input:

Output:

Landscape Position Continued– Median and ImageMagick

Highlighting ridges with 250ft buffer (on 2.5ft DEM) with just ImageMagick:


convert lscape_posit.png -median 100 median100.png
composite -compose difference lscape_posit.png median100.png difference_median100.png


Input:

Output:

BTW, median calculations of this size are slow, even in ImageMagick.

Landscape Position Continued– absolutely relative position calculation

I apologize in advance– this first post will be heavy on code, short on explanation. Landscape position, e.g. previous posts, can be trivial to calculate, but to make the calculations scalable to a large area, some batching is necessary. In this case, instead of a McNab index, we’re calculating the traditional GIS landscape position. Enter my favorite non-geographic tool, PovRay… . To the difference between, we’ll use PovRay in combination with imagemagick:

#include "transforms.inc"
#version 3.6;

#declare widthx=3425;
#declare heighty=1707.5;

#declare Aperture=300;
#declare Laps=20;

#declare Ind=1-pow(1-clock,1.2);
#declare PosX=cos(2*pi*Laps*Ind)*Aperture*Ind;
#declare PosY=sin(2*pi*Laps*Ind)*Aperture*Ind;

#declare Camera_Location = ;
#declare Camera_Lookat = ;

camera {
	orthographic
	location Camera_Location
	look_at Camera_Lookat
	right widthx*x
	up heighty*y
}

background {color  }

union {

	height_field {
		png "dem.png"
		scale ;
		translate ;
	}

	pigment {
		gradient x color_map {
			[0 color rgb 1]
			[1 color rgb 0]
			}

		scale 
		Reorient_Trans(x, Camera_Lookat-Camera_Location)
		translate Camera_Location
	}

	finish {ambient 1}
}

povray +Ilscape_posit.pov +Olscape_posit.png +FN16 +W1370 +H683 +KFI1 +KFF99 -D
convert lscape_posit??.png -average average.png
composite -compose difference lscape_posit.png average.png difference.png

In truth, I think I could do all of this in imagemagick, but it might not be fast enough. More testing to follow... .