Running statistics using RunningStats (how to run on Python using SWIG and Distutils)

There is a more accurate method for calculating focal statistics than usually implemented in numerical packages. This is documented here and the C++ code is download from here

The SWIG setup file is simple:



which you run using:



The Distutils setup script is:



compile using "python setup.py install"

Then in python, a simple implementation of calculating the global stats from an array would be something like:



 To show the convergence of locals means on the global mean, you might use the RunningStats class in the following way:


 Finally, to compute running averages on the array you might use something like the following, where the function to window the data efficiently was downloaded from here


Obviously, as well as the mean, the above could be used for variance, standard deviation, kurtosis and skewness.

Pymix fix: short note on compiling PyMix on Ubuntu

A fix for users of Python 2.6 or 2.7 when installing PyMix:

How to pass arguments to python programs ...

... the right way (in my opinion). This is the basic template I always use. Not complicated, but the general formula might be useful to someone.

It allows you to pass arguments to a program with the use of flags, so they can be provided in any order. It has a 'help' flag (-h), will raise an error if no arguments are passed.

In this example one of the inputs (-a) is supposed to be a file. If the file is not provided, it uses Tkinter to open a simple 'get file' window dialog.

It then parses the other arguments out. Here I have assumed one argument is a string, one is a float, and one is an integer. If any arguments are not provided, it creates variables based on default variables. It prints everything it's doing to the screen.

Simple Python script to write a list of coordinates to a .geojson file

This has probably been done better elsewhere, but here's a really simple way to write a geojson file programmatically to display a set of coordinates, without having to bother learning the syntax of yet another python library. Geojson is a cool little format, open source and really easy.

Below, arrays x and y are decimal longitude and latitude respectively. The code would have to be modified to include a description for each coordinate, which could be included easily on line 15 if in the form of an iterable list.

Doing Spectral Analysis on Images with GMT

Mostly an excuse to play with, and show off the capabilities of, the rather fabulous Ipython notebook! I have created an Ipython notebook to demonstrate how to use GMT to do spectral analysis on images.

This is demonstrated by pulling wavelengths from images of rippled sandy sediments.

GMT is command-line mapping software for Unix/Linux systems. It's not really designed for spectral analysis in mind, especially not on images, but it is actually a really fast and efficient way of doing it! 

This also requires the convert facility provided by Image Magick. Another bodacious tool I couldn't function without.

My notebook can be found here.

Fantastic images of ripples provided by http://skeptic.smugmug.com

Destination point given distance and bearing from start point

You have starting points:

xcoordinate, a list of x-coordinates (longitudes)

ycoordinate, a list of y-coordinates (latitudes)

num_samples, the number of samples in the plane towards the destination point

bearings:

heading, a list of headings (degrees)

and distances:

range, a list of distances in two directions (metres)

This is how you find the coordinates (x, y) along the plane defining a starting point and two end points, in MATLAB.

Coordinate system conversion using cs2cs from PROJ

I hate geodesy. It's confusing, multifarious, and ever-changing (and I call myself a scientist!). But we all need to know where we, and more importantly our data, are. I have need of a single tool which painlessly does all my conversions, preferably without me having to think too much about it.

You could use an online tool, but when you're scripting (and I'm ALWAYS scripting) you want something a little less manual. Enter cs2cs from the PROJ.4 initiative. It will convert pretty much anything to anything else. Perfect.

In Ubuntu:
sudo apt-get install proj-bin

In Fedora I searched yum for proj-4

If the man page is a little too much too soon to take in, this wonderful page helps you decide which flags and parameters to use, with a rudimentary understanding of what you want

For example, if I want to convert coordinates in Arizona Central State Plane to WGS84 Latitude/Longitude, I use this:



(the State Plane Coordinate System is new to me, but my view of the whole field of geodesy is 'let's make things even more complicated, sorry, accurate!')

Let's decompose that a bit:



tells it I want decimal degrees with 6 values after the decimal point (you could leave this out if you want deg, min, sec)



says use the Transverse Mercator projection



are all parameters related to the input coordinates in Arizona Central State Plane (the 'to_meter' is a conversion from feet, the default unit, to metres)



are the parameters related to the output coordinates in Latitude/Longitude

infile is a two column list of points e.g.
2.19644131000e+005 6.11961823000e+005
2.17676234764e+005 6.11243565478e+005
2.19763457634e+005 6.11234534908e+005
2.19786524782e+005 6.11923555789e+005
2.19762476867e+005 6.11378246389e+005

outfile will be created with the results:
-111.846501    36.517758 0.000000
-111.868478    36.511296 0.000000
-111.845175    36.511202 0.000000
-111.844912    36.517412 0.000000
-111.845185    36.512498 0.000000

By default it will always give you the third dimension (height) even if you didn't ask for it, like above, so I tend to trim it by asking awk to give me just the first 2 columns separated by a tab:



-111.846501      36.517758
-111.868478      36.511296
-111.845175      36.511202
-111.844912      36.517412
-111.845185      36.512498

Compiling points2grid on Fedora with g++

Gridding large [x,y,z] point cloud datasets. Matlab and python have both failed me with their inefficient use of memory. So I'm having to turn to other tools.

Points2Grid is a tool, written in C++, which promises to do the job in a simple intuitive way. It accepts ascii and las (lidar) formats and outputs gridded data, plus gridded data densities.

First install LIBLAS which is Lidar data translation libraries. The recommended version is 1.2.1. I installed liblas 1.6.1 because it was the oldest version on the site that had a working link. I installed this by issuing:

Make a note of where the libraries are installed. On my system: /usr/local/include/liblas/

Download points2grid from here. 

In Interpolation.h, adjust MEM_LIMIT variable in the  file to reflect the memory limitations of your system. From the README guide, 

"As a rule of thumb, set the MEM_LIMIT to be (available memory in bytes)/55. For jobs that generate

grids cells greater than the MEM_LIMIT, points2grid goes "out-of-core", resulting in significantly worse performance."

Find this value using:

I then had to point the code to where the liblas directories had been installed. In Interpolation.cpp, replace:

It should then install with a simple 'make' command. Happy gridding!

Automatic Clustering of Geo-tagged Images. Part 2: Other Feature Metrics

In the last post I took a look at a simple method for trying to automatically cluster a set of images (of Horseshoe Bend, Glen Canyon, Arizona). Some of those picture were of the river (from different perspectives) and some of other things nearby.

With a goal of trying to separate the two classes, the results were reasonably satisfactory with image histograms as a feature metric, and a little user input. However, it'd be nice to get a more automated/objective way to separate the images into two discrete clusters. Perhaps the feature detection method requires a little more scrutiny?

In this post I look at 5 methods using the functionality of the excellent Mahotas toolbox for python downloaded here:

1) Image moments (I compute the first 5)
2) Haralick texture descriptors which are based on the co-occurrence matrix of the image
3) Zernike Moments
4) Threshold adjacency statistics (TAS)
5) Parameter-free threshold adjacency statistics (PFTAS)

The code is identical to the last post except for the feature extraction loop. What I'm looking for is as many images of the river at either low distances (extreme left in the dendrogram) or very high distances (extreme right in the dendrograms). These clusters have been highlighted below. There are 18 images of the river bend.

In order of success (low to high):

1) Zernike (using a radius of 50 pixels and 8 degrees, but other combinations tried):

 This does a poor job, with no clusters at either end of the dendrogram.

2) Haralick:

There is 1 cluster of 4 images at the start

3) Moments:

There is 1 cluster of 6 images at the start, and 1 cluster of 2 at the end

4) TAS:

There is 1 cluster of 7 images at the start, and 1 cluster of 5 near the end, and both clusters are on the same stem

5) PFTAS:

The cluster at the end contains 16 out of 18 images of the river bend at the end. Plus there are no user-defined parameters. My kind of algorithm. The winner by far!

So there you have it, progress so far! No methods tested so far is perfect. There may or may not be a 'part 3' depending on my progress, or lack of!

Automatic Clustering of Geo-tagged Images. Part 1: using multi-dimensional histograms

There's a lot of geo-tagged images on the web. Sometimes the image coordinate is slightly wrong, or the scene isn't quite what you expect. It's therefore useful to have a way to automatically download images from a given place (specified by a coordinate) from the web, and automatically classify them according to content/scene.

In this post I describe a pythonic way to:
1) automatically download images based on input coordinate (lat, long)
2) extract a set features from each image
3) classify each image into groups
4) display the results as a dendrogram

This first example, 2) is achieved using a very simple means, namely the image histogram of image values. This doesn't take into account any texture similarities or connected components, etc. Nonetheless, it does a reasonably good job at classifying the images into a number of connected groups, as we shall see.

In subsequent posts I'm going to play with different ways to cluster images based on similarity, so watch this space!

User inputs: give a lat and long of a location, a path where you want the geo-tagged photos, and the number of images to download




First, import the libraries you'll need. Then interrogates the website and downloads the images.



As you can see the resulting images are a mixed bag. There's images of the river bend, the road, the desert, Lake Powell and other random stuff. My task is to automatically classify these images so it's easier to pull out just the images of the river bend

The following compiles a list of these images. The last part is to sort which is not necessary but doing so converts the list into a numpy array which is.  Then clustering of the images is achieved using Jan Erik Solem's rather wonderful book 'Programming Computer Vision with Python'. The examples from which can be downloaded here. Download this one, then this bit of code does the clustering:

The approach taken is to use hierarchical clustering using a simple euclidean distance function. This bit of code does the dendogram plot of the images:

which in my case looks like this (click on the image to see the full extent):

It's a little small, but you can just about see (if you scroll to the right) that it does a good job at clustering images of the river bend which look similar. The single-clustered, or really un-clustered, images to the left are those of the rim, road, side walls, etc which don't look anything like the river bend. 

Next, reorder the image list in terms of similarity distance (which increases in both the x and y directions of the dendrogram above)

Which gives me:

As you can see they are all images of the river bend, except the 2nd from the left on the top row, which is a picture of a shrub. Interestingly, the pattern of a circular shrub surrounded by a strip of sand is visually similar to the horse shoe bend!!

However, we don't want to include it with images of the river, which is why a more sophisticated method that image histogram is required to classify and group similar image ... the subject of a later post.

Alpha Shapes in Python

Alpha shapes include convex and concave hulls. Convex hull algorithms are ten a penny, so what we're really interested in here in the concave hull of an irregularly or otherwise non-convex shaped 2d point cloud, which by all accounts is more difficult.

The function is here:



The above is essentially the same wrapper as posted here, except a different way of reading in the data, and providing the option to specify a probe radius. It uses Ken Clarkson's C-code, instructions on how to compile are here

An example implementation:



The above data generation is translated from the example in a matlab alpha shape function . Then:



Which produces (points are blue dots, alpha shape is red line):

and on a larger data set (with radius=1, this is essentially the convex hull):

and a true concave hull:

Installing Mb-System on Fedora 17

MB-System is the only fully open-source software I am aware of for swath sonar data including multibeam and sidescan backscatter. As such, I'm really hoping using it works out for me!

Compiling it on a fresh install of Fedora 17 wasn't that straight-forward (took me about 5 hours to figure all this out), so I thought I'd put together a little how-to for the benefit of myself, and hopefully you.

I've broken it down into a few stages so some sort of semblance of order can be imposed on this random collection of installation notes.

Step 1: download the software from ftp server here

In the terminal:

For me this created a folder called 'mbsystem-5.3.1982'

Step 2: install pre-requisites

A) Generic Mapping Tools (GMT). You need to go to the Downloads tab, then click on the INSTALL FORM link. This will take you to a page where you can download the installation script, and fill in an online form of your parameter settings which you can submit and save the resulting page as a textfile which becomes the input to the bash script. Sounds confusing, but it's not - the instructions on the page are adequate. I opted to allow GMT to install netCDF for me. Then in the terminal I did this:

B) Perl parallel forkmanager:

C) X11

D) openmotif

E) fftw

F) ghostview - I had to install this indirectly using kde3:

(Note to Fedora - why oh why oh why are B, C, and F above not installed by default!?)

G)  OTPSnc tidal prediction software: download from here

untar, and cd to the directory

first double check that ncdump and ncgen are installed (which ncdump ncgen)

then edit the makefile so it reads:

then in the terminal issue:

Hopefully this compiles without errors, then I moved them to a executable directory:

H) openGL libraries (again I had to do this indirectly, this time through mesa):

Step 3: prepare mbsystem makefiles

cd mbsystem-5.3.1982/

You have to go in and point install_makefiles to where all your libraries are. This is time-consuming and involves a lot of ls, which, and whereis!

Here's a copy of the lines I edited in my install_makefiles parameters:

Then in the terminal

Step 3: install mbsystem

first I had to do this (but you may not need to)

I then updated my ~/.bashrc so the computer can find all these lovely new files:

Test the install by issuing



you should see the manual pages come up. Right, now down to the business of multibeam analysis!

Alternative to bwconvhull

Matlab's bwconvhull is a fantastic new addition to the image processing library. It is a function which, in its full capacity, will return either the convex hull of the entire suite of objects in a binary image, or alternatively the convex hull of each individual object. I use it, in the latter capacity, for filling multiple large holes (holes large enough that imfill has no effect) in binary images containing segmented objects 

For those who do not have a newer Matlab version with bwconvhull included, I have written this short function designed to be an alternative to the usage

 P= bwconvhull(BW,'objects');

and here it is:



It uses the IP toolbox function regionprops, around for some time, to find the convex hull images of objects present, then inscribes them onto the output image using the bounding box coordinates. Simples!

Ellipse Fitting to 2D points, part 3: Matlab

The last of 3 posts presenting algorithms for the ellipsoid method of Khachiyan for 2D point clouds. This version in Matlab will return: 1) the ellipse radius in x, 2) the ellipse radius in y, 3) the centre in x, 4) the centre in y, and 5) orientation of the ellipse, and coordinates of the ellipse

Ellipse Fitting to 2D points, part 1: R

The first of 3 posts presenting algorithms for the ellipsoid method of Khachiyan for 2D point clouds. This version in R will return: 1) the ellipse radius in x, 2) the ellipse radius in y, 3) the centre in x, 4) the centre in y, and 5) orientation of the ellipse, and coordinates of the ellipse

Non-linear least squares fit forced through a point

Sometimes you want a least-squares polynomial fit to data but to constrain it so it passes through a certain point. Well ... Inputs: x, y (data); x0 and y0 are the coordinates of the point you want it to pass through, and n is the degree of polynomial

Program for spatial correlogram of an image, in C

This program will read an image, convert it to greyscale, crop it to a square, and calculate it's spatial (1D) autocorrelation function (correlogram), printing it to screen.

It uses OpenCV image processing libraries - see here and here

Compilation instructions (in Ubuntu):

 sudo apt-get install libcv-dev libcv1 libcvaux-dev libcvaux1 libhighgui-dev libhighgui1  
 pkg-config --libs opencv; pkg-config --cflags opencv  
 PKG_CONFIG_PATH=/usr/lib/pkgconfig/:${PKG_CONFIG_PATH}  
 export PKG_CONFIG_PATH  
 gcc `pkg-config --cflags --libs opencv` -o Rspat Rspat.c  

Usage in the terminal:

 ./Rspat IMG_0208.JPG   

The program:

Updated Power Spectrum from an Image

A while back I posted on using ImageMagick to do a full frequency-decomposition on an image, and create a power spectrum. It had several inaccuracies. Here's an update with how to do it better

First, you need to install a version of IM which is HDRI enabled (see here)

# test for valid version of Image Magick

 im_version=`convert -list configure | \  
  sed '/^LIB_VERSION_NUMBER /!d; s//,/; s/,/,0/g; s/,0*\([0-9][0-9]\)/\1/g'`  
 [ "$im_version" -lt "06050407" ] && errMsg "--- REQUIRES IM VERSION 6.5.4-7 OR HIGHER ---"  

# test for hdri enabled

 hdri_on=`convert -list configure | grep "enable-hdri"`  
 [ "$hdri_on" = "" ] && errMsg "--- REQUIRES HDRI ENABLED IN IM COMPILE ---"  

#read $file, load and converting to greyscale

 convert $file -colorspace gray in_grey.tiff  

# crop image to smallest dimension

 min=`convert in_grey.tiff -format "%[fx: min(w,h)]" info:`  

# crop image, and convert to png (why?)

 convert in_grey.tiff -background white -gravity center -extent ${min}x${min} in_grey_crop.png  

# get data, find mean as proportion, convert mean to percent, and de-mean image

 data=`convert in_grey_crop.png -verbose info:`  
 mean=`echo "$data" | sed -n '/^.*[Mm]ean:.*[(]\([0-9.]*\).*$/{ s//\1/; p; q; }'`  
 m=`echo "scale=1; $mean *100 / 1" | bc`  
 convert in_grey_crop.png -evaluate subtract ${m}% in_grey_crop_dm.png  

# pad to power of 2

 p2=`convert in_grey_crop_dm.png -format "%[fx:2^(ceil(log(max(w,h))/log(2)))]" info:`  
 convert in_grey_crop_dm.png -background white -gravity center -extent ${p2}x${p2} in_grey_crop_dm_pad.png  

# fft

 convert in_grey_crop_dm_pad.png -fft in_ft.png  

# get spectrum

 convert in_ft-0.png -contrast-stretch 0 -evaluate log 10000 in_spec.png  

... but it's still not quite right. Any tips?

Installing NCVIEW for Ubuntu

Recently I ran into problems installing NCVIEW, which is a tool for efficiently viewing the contents of a NetCDF data file. The problem was something to do with the X terminal. So here's how to install it from source:

 wget ftp://cirrus.ucsd.edu/pub/ncview/ncview-1.93g.tar.gz  
 sudo apt-get install libxaw7-dev  
 tar -xvvf ncview-1.93g.tar.gz  
 cd ncview-1.93g/  
 ./configure --x-libraries=/usr/X11R6/lib  
 make  
 make install