How to call PostGIS ST_* functions in GeoDjango

In short: use django queryset extra()

I needed to build a web interface using GeoDjango for querying an existing (huge) PostGIS database table which has a geometry column:

the_geom = models.PointField(srid=4326)

python manage.py inspectdb conveniently created models.py for my django application with the correct class Meta db_table set to the legacy database table name. I just added ‘managed = False‘ so that django does not try to manage the table in any way.

Basically, my application requires the user to upload a KML file containing geometries of interest (e.g. POINT, LINESTRING, etc.), input a distance value in metres, and optional date range parameters.

The application then returns all the points from the database table which are within the user-specified distance and geometries of interest in KML format.

GeoDjango provides two spatial lookup APIs for distance calculation, dwithin and distance. dwithin runs many times faster than distance but it requires the distance to be given in degrees.

PostGIS documentation suggested using geography() and ST_DWithin() for best performance. But, how to I call postGIS geography() in a where-clause from my django application?

I have considered writing a custom manager and using raw() but raw() returns a RawQueryset, which unlike the ordinary Queryset, does not have filter() method (thus chaining of results is not possible). Besides, I would like to make use of postGIS ST_AsKml() on the query result.

Putting everything together:

MyModel.objects.extra(
    select={'the_geom_wkt': 'ST_AsText(the_geom)'},
    where=[
        "ST_DWithin(geography(the_geom), ST_GeographyFromText('SRID=4326;%s'),%%s)"%%geometry_from_kml
    ],
    params=[distance_in_metres],
)

How to make a circle geometry in KML or WKT

There is no circle geometry in KML.

In remote sensing, some algorithms e.g. MODIS wild fire detection, etc. produce sparse point data. But, I find it more meaningful to visualise these point data on a map, Google Earth layer for instance, as circles whose radii are based on the satellite imagery resolution.

As suggested in the KML FAQ, I approximate the circle geometry using regular n-sided polygon. The Python code below is originates from las3rjock’s posting on StackOverflow. I have adapted it slightly to accept input in metres, and to produce KML fragment and WKT. It is good enough when radius is small compared to the radius of the Earth.

Alternatively, get KML circle which can also do stars, or use the online generator.

from math import asin,cos,pi,sin
import re
class GeomHelper:
    rEarth_km = 6371.01 # Earth's average radius in km
    epsilon = 0.000001 # threshold for floating-point equality
    @staticmethod
    def deg2rad(angle):
        return angle*pi/180
    @staticmethod
    def rad2deg(angle):
        return angle*180/pi
    @staticmethod
    def pointRadialDistance(lat1_deg, lon1_deg, bearing_deg, distance_m):
        """
        Return final coordinates (lat2,lon2) [in degrees] given initial coordinates
        (lat1,lon1) [in degrees] and a bearing [in degrees] and distance [in m]
        """
        rlat1 = GeomHelper.deg2rad(lat1_deg)
        rlon1 = GeomHelper.deg2rad(lon1_deg)
        rbearing = GeomHelper.deg2rad(bearing_deg)
        rdistance = distance_m / (GeomHelper.rEarth_km*1000) # normalize linear distance to radian angle
        rlat = asin( sin(rlat1) * cos(rdistance) + cos(rlat1) * sin(rdistance) * cos(rbearing) )
        if cos(rlat) == 0 or abs(cos(rlat)) < GeomHelper.epsilon: # Endpoint a pole
            rlon=rlon1
        else:
            rlon = ( (rlon1 - asin( sin(rbearing)* sin(rdistance) / cos(rlat) ) + pi ) % (2*pi) ) - pi
        lat = GeomHelper.rad2deg(rlat)
        lon = GeomHelper.rad2deg(rlon)
        return (lat, lon)

    @staticmethod
    def make_polygon(lon_deg, lat_deg, radius_m, sides=8, format='wkt'):
        theta = 360/sides
        if format is 'wkt':
            separator = ' '
        else:
            separator = ','
        result = []
        for i in range(sides):
            lat, lon = GeomHelper.pointRadialDistance(lat_deg, lon_deg, i*theta, radius_m)
            result.append("%f%s%f" % (lon, separator, lat))
        # Close the polygon
        result.append(result[0])
        if format is 'wkt':
            return 'POLYGON((%s))' % ','.join(result)
        else:
            # KML
            return '%s' % ' '.join(result)
    @staticmethod
    def make_polygon_from_wkt(point, radius_m, **kwargs):
        """ 'point' is a string like POINT(30,-25.45) """
        m = re.search('^point\((\S+) ([^)]+)\)$', point, re.I)
        if m:
            return GeomHelper.make_polygon(float(m.group(1)), float(m.group(2)), radius_m, **kwargs)
        else:
            raise ValueError("Invalid WKT POINT: %s" % point)
    @staticmethod
    def make_polygon_from_kml(point, radius_m, **kwargs):
        """ 'point' is a string like <POINT><coordinates>27.113,-25.389</coordinates></POINT> """
        m = re.search('<point><coordinates>([^,]+),([^<]+)', point, re.I)
        if m:
            return GeomHelper.make_polygon(float(m.group(1)), float(m.group(2)), radius_m, **kwargs)
        else:
            raise ValueError("Invalid KML POINT: %s" % point)

GeoTIFF and Python GDAL

If you program in Python, then the Python binding to the Geospatial Data Abstraction library (GDAL) and Numpy are indispensable tools for reading, writing and manipulating raster image in GeoTIFF format.

This summarizes my brief experience working with them.

A GeoTIFF image file may contain one or more bands. Each band can be viewed as a separate 2D data array but every band in the same file will have the same array dimension (or shape, in Numpy lingo).

All image pixels in a band (data element of array) are of the same data type, e.g. byte, int16, etc., but one band may be of different data type from another in the same file.

By convention, if a GeoTIFF file contains a RGB raster image, it will have 3 bands (4 if it is RGBA, with the addition of alpha channel for pixel transparency) of byte data — band 1 for red, 2 for green and 3 for blue.

GDAL API

To read GeoTIFF metadata (georeference information)

from osgeo import gdal
ds = gdal.Open("sample.tif", gdal.GA_ReadOnly)
(X, deltaX, rotation, Y, rotation, deltaY) = ds.GetGeoTransform()
srs_wkt = ds.GetProjection()
Nx = ds.RasterXSize
Ny = ds.RasterYSize
ds = None

GeoTransform() returns a tuple where (X, Y) are corner coordinates of the image indicating the origin of the array, i.e. data element[0,0]. But, which corner?

If deltaX is positive, X is West. Otherwise, X is East. If deltaY is positive, Y is South. Otherwise, Y is North. In other words, when both deltaX and deltaY is positive, (X, Y) is the lower-left corner of the image.

It is also common to have positive deltaX but negative deltaY which indicates that (X, Y) is the top-left corner of the image.

There are several standard notations for SRS, e.g. WKT, Proj4, etc. GetProjection() returns a WKT string. The meaning and unit-of-measurement for  X, Y, deltaX and deltaY are based on the spatial reference system (SRS) of the image. For example, if the SRS is Stereographic projection, they are in kilometres.

To convert from WKT notation to Proj4

from osgeo import osr
srs = osr.SpatialReference()
srs.ImportFromWKT(srs_wkt)
srs_proj4 = srs.ExportToProj4()

To load a multi-band GeoTIFF image file as Numpy arrays

from osgeo import gdal
arys = []
ds = gdal.Open('sample.tif', gdal.GA_ReadOnly)
for i in xrange(1, ds.RasterCount+1):
arys.append(ds.GetRasterBand(i).ReadAsArray())
# for a single-band GeoTIFF
ary = ds.GetRasterBand(1).ReadAsArray()

Unlike mathematical convention for coordinates (X, Y), Numpy array indexing uses [y-offset][x-offset].

To save a 2-D Numpy array as 1-band GeoTIFF file

from osgeo import gdal
 # data exists in 'ary' with values range 0 - 255
 # Uncomment next line if ary[0][0] is upper-left corner
 #ary = numpy.flipup(ary)
 Ny, Nx = ary.shape
 driver = gdal.GetDriverByName("GTiff")
 ds = driver.Create('output.tif', Ny, Nx, 1, gdal.GDT_Byte)
 ds.SetGeoTransform( ... ) # define GeoTransform tuple
 ds.GetRasterBand(1).WriteArray(ary)
 ds = None

To slice from a bigger image (array)

subary = ary[4:6, 5:8]

To change coordinates from one SRS to another

Use pyproj
For example, to convert 2 lists of longitudes and latitudes to X, Y coordinates in polar stereographic projection:

import pyproj
# Specify target SRS
pp = pyproj.Proj('+proj=stere +lat_ts=-60.0 +lat_0=-90.0 +lon_0=28.0 +k_0=1.0')

# Origin of 2D-data array: longitudes[0,0] latitudes[0,0]
X, Y = pp(longitudes, latitudes)

# To reverse the operation
longitudes, latitudes = pp(X, Y, inverse=True)

How to follow a big list of Twitter users

Say you have a big list of Twitter users you want to follow (example), you can add the Twitter username one-by-one or you can leverage the power of command line on Linux. Here’s how:

1. Download the twidge binary appropriate for your Linux
2. Rename the downloaded binary file as ‘twidge’ and make it executable i.e. chmod u+x twidge
3. twidge setup    # Follow instructions given to authorize twidge application to access your Twitter account
4. Using the example above, copy and paste the list of Twitter usernames and save the file as ‘userlist’
5. cat userlist | sed ‘s/@/\n/g’ | while read U X; do if [ -z "$U" ]; then continue; fi; twidge follow $U; done

How to parse HTML using Python HTMLParser

A sample code snippet on how to use the Python module HTMLParser to extract a well-formed HTML document for multiple

<input name="fileIDs" value="123456" />

Code snippet:

import urllib, urllib2

from HTMLParser import HTMLParser

class MyHTMLParser(HTMLParser):

    def __init__(self, fh):

        """

        {fh} must be an input stream returned by open() or urllib2.urlopen()

        """

        HTMLParser.__init__(self)

        self.fileids = []

        self.feed(fh.read())

    def handle_starttag(self, tag, attrs):

        if tag == 'input':

            attrD = dict(attrs)

            if attrD['name'] == 'fileIDs':

                self.fileids.append(attrD['value'])

    def get_fileids(self):

        return self.fileids

opener = urllib2.build_opener(urllib2.HTTPHandler(debuglevel=1))

opener.addheaders = [('User-agent', 'Mozilla/5.0')]

response = opener.open("http://www.example.com/200.html")

myparser = MyHTMLParser(response)

Reference:

• http://www.hellboundhackers.org/articles/841-using-python-39;s-htmlparser-class.html

How to migrate 32-bit Cyrus IMAPD mailboxes to 64-bit

Thanks to Viktor Petersson’s blog for the crucial hints. I adapted his procedure for FreeBSD slightly because I was moving the cyrus-imapd service from an old 32-bit SuSE Linux Enterprise Server version 9 (SLES9) to a new 64-bit server running CentOS 5.5.

On old server,

rsync -aP /var/lib/imap root@mynewserver:/var/lib/imap
rsync -aP /var/spool/imap root@mynewserver:/var/spool/imap

Stop Cyrus service on both old and new servers.

rsync -aP --delete /var/lib/imap root@mynewserver:/var/lib/imap
rsync -aP --delete /var/spool/imap root@mynewserver:/var/spool/imap
sudo -u cyrus /usr/local/cyrus/bin/ctl_mboxlist -d > ~/mboxlist.txt

Copy mboxlist.txt to new server.

On the new server,

sudo -u cyrus rm /var/lib/imap/db/* \
/var/lib/imap/db.backup1/* \
/var/lib/imap/db.backup2/* \
/var/lib/imap/deliver.db \
/var/lib/imap/tls_sessions.db \
/var/lib/imap/mailboxes.db
sudo -u cyrus /usr/lib/cyrus-imapd/ctl_mboxlist -u < ~/mboxlist.txt
sudo -u cyrus /usr/lib/cyrus-imapd/ctl_cyrusdb -r
sudo -u cyrus /usr/lib/cyrus-imapd/tls_prune
sudo -u cyrus /usr/lib/cyrus-imapd/ctl_cyrusdb -c
sudo -u cyrus /usr/lib/cyrus-imapd/cyr_expire -E 3
sudo -u cyrus mkdir /var/lib/imap/rpm/
cd /var/spool/imap/users
for U in * ; do 
  sudo -u cyrus /usr/lib/cyrus-imapd/reconstruct -r -f "user.$U"; 
done

Start the service on the new server.

Six stages of a project

Applicable to any project involving more than 5 people. Author unknown.

1. Enthusiasm
2. Disillusionment
3. Panic
4. Persecution of the innocent
5. Praise of the bystander

A collection of Latin quotes I like

I think I had my first encounter with Latin quotes/proverbs from reading Asterix comics (thanks Matthias). Perhaps it’s simply because I don’t know Latin that I like them being so concise and imperative.

Ponte facto Caesar transit.

“The ablative absolute,” his Latin teacher had explained, “is used by men of action who don’t want to waste words. Ponte facto Caesar transit. The bridge built, Caesar crossed it. Consider how effective this is. No bothering with who built the bridge or at what cost. The bridge was built, as a bridge should be, and Caesar crossed it.” (source)

Nullius in verba

“Take nobody’s word for it” (source)

Acta est fabula

It’s all over (lit. the drama has been acted out) (source)

Alea jacta est

The die is cast

Audaces fortuna juvat

Fortune favors the bold

Auri sacra fames

The cursed hunger for gold

Bis repetita placent

The things that pleases are repeated again and again

Cognito ergo sum

I think therefore I am

Ab igne ignem capere

“To light a fire with a fire.” (source)

How to install Novell Groupwise 8 on 32-bit Ubuntu 9.10, 10.4, 10.10

The following Ubuntu packages are required:

1. alien
2. sun-java6-bin
3. sun-java6-jre

Next,

1. Download the latest ZIP archive from http://gwclient.provo.novell.com/client/GW801LinuxClient.zip
2. Uncompress the ZIP archive to extract the following two RPM files. The actual  RPM version may vary depending on the Groupwise ZIP package downloaded.
   1. novell-groupwise-client-8.0.1-88138.i586.rpm
   2. novell-groupwise-gwcheck-8.0.1-88138.i586.rpm
3. sudo alien novell-groupwise-client-8.0.1-88138.i586.rpm to produce novell-groupwise-client_8.0.1-88139_i386.deb
4. sudo dpkg -i novell-groupwise-client_8.0.1-88139_i386.deb
5. Edit /usr/bin/groupwise to replace a line near the beginning that goes, “export LD_LIBRARY_PATH=…” with the following:

# JRE_HOME may vary depending on the exact version you have installed
JVMHOME=$(ls -d /usr/lib/jvm/java-6-sun-* | tail -1)
if [ -z "$JVMHOME" ]; then
 echo "Install sun-java6-jre and sun-java6-bin packages first" >&2
 exit 1
fi
export JRE_HOME=$JVMHOME/jre
export LD_LIBRARY_PATH=$JRE_HOME/lib/i386:$JRE_HOME/lib/i386/client:$LD_LIBRARY_PATH

 

24 Timeless quotes on politics

From a forwarded email:

The adage that “ the only thing Man learns from History is that Man learns nothing from History” seems to be backed up by sayings of famous people, going as far back as 430 BC!

1. In my many years I have come to a conclusion that one useless man is a shame, two is a law firm and three or more is a congress. – John Adams
2. If you don’t read the newspaper you are uninformed, if you do read the newspaper you are misinformed. – Mark Twain
3. Suppose you were an idiot. And suppose you were a member of Congress. But then I repeat myself. — Mark Twain
4. I contend that for a nation to try to tax itself into prosperity is like a man standing in a bucket and trying to lift himself up by the handle … – Winston Churchill
5. A government which robs Peter to pay Paul can always depend on the support of Paul. – George Bernard Shaw
6. A liberal is someone who feels a great debt to his fellow man, which debt he proposes to pay off with your money. – G. Gordon Liddy
7. Democracy must be something more than two wolves and a sheep voting on what to have for dinner. – James Bovard, Civil Libertarian (1994)
8. Foreign aid might be defined as a transfer of money from poor people in rich countries to rich people in poor countries. – Douglas Casey, Classmate of Bill Clinton at Georgetown University
9. Giving money and power to government is like giving whiskey and car keys to teenage boys. – P.J. O’Rourke, Civil Libertarian
10. Government is the great fiction, through which everybody endeavors to live at the expense of everybody else. – Frederic Bastiat, French Economist (1801-1850)
11. Government’s view of the economy could be summed up in a few short phrases: If it moves, tax it. If it keeps moving, regulate it. And if it stops moving, subsidize it. – Ronald Reagan (1986)
12. I don’t make jokes. I just watch the government and report the facts. – Will Rogers
13. If you think health care is expensive now, wait until you see what it costs when it’s free! – P.J. O’Rourke
14. In general, the art of government consists of taking as much money as possible from one party of the citizens to give to the other. – Voltaire (1764)
15. Just because you do not take an interest in politics doesn’t mean politics won’t take an interest in you! – Pericles (430 B.C.)
16. No man’s life, liberty, or property is safe while the legislature is in session. – Mark Twain (1866)
17. Talk is cheap…except when Congress does it. – Anonymous
18. The government is like a baby’s alimentary canal, with a happy appetite at one end and no responsibility at the other. – Ronald Reagan
19. The inherent vice of capitalism is the unequal sharing of the blessings. The inherent blessing of socialism is the equal sharing of misery. – Winston Churchill
20. The only difference between a tax man and a taxidermist is that the taxidermist leaves the skin. – Mark Twain
21. The ultimate result of shielding men from the effects of folly is to fill the world with fools. – Herbert Spencer, English Philosopher (1820-1903)
22. There is no distinctly native American criminal class…save Congress. – Mark Twain
23. What this country needs are more unemployed politicians. – Edward Langley, Artist (1928-1995)
24. A government big enough to give you everything you want, is strong enough to take everything you have. – Thomas Jefferson