(sorry, I’ve been busy in the last few months, hence no blogposts)

Update: almost at the same time I originally posted this, Sean set me up on the http://gispython.org blog, so I’ve moved this post there.

PyWPS is a neat Python package supporting the OGC Web Processing Service standard.  Basic setup and configuration can be found in the documentation, or Tim’s useful post.

I’ve been working on a demo to expose the OGR Python bindings for geoprocessing (buffer, centroid, etc.).

Here’s an example process to buffer a geometry (input as WKT), and output either GML, JSON, or WKT:

from pywps.Process import WPSProcess
import osgeo.ogr as ogr

class Buffer(WPSProcess):
 def __init__(self):
  WPSProcess.__init__(self,
  identifier='buffer',
  title='Buffer generator',
  metadata=['http://www.kralidis.ca/'],
  profile='OGR / GEOS geoprocessing',
  abstract='Buffer generator',
  version='0.0.1',
  storeSupported='true',
  statusSupported='true')

  self.wkt = self.addLiteralInput(identifier='wkt', \
   title='Well Known Text', type=type('string'))
  self.format = self.addLiteralInput(identifier='format', \
   title='Output format', type=type('string'))
  self.buffer = self.addLiteralInput(identifier='buffer', \
   title='Buffer Value', type=type(1))
  self.out = self.addLiteralOutput(identifier='output', \
   title='Buffered Feature', type=type('string'))

 def execute(self):
  buffer = ogr.CreateGeometryFromWkt( \
   self.wkt.getValue()).Buffer(self.buffer.getValue())
  self.out.setValue(_genOutputFormat(buffer, self.format.getValue()))
  buffer.Destroy()

 def _setNamespace(xml, prefix, uri):

  return xml.replace('>', ' xmlns:%s="%s">' % (prefix, uri), 1)

 def _genOutputFormat(geom, format):

 if format == 'gml':

  return _setNamespace(geom.ExportToGML(), 'gml', \
     'http://www.opengis.net/gml')

 if format == 'json':

  return geom.ExportToJson()

 if format == 'wkt':

  return geom.ExportToWkt()

Notes:

• _setNamespace is a workaround, as OGR’s ExportToGML doesn’t declare a namespace prefix / uri in the output, which would make the ExecuteResponse XML choke parsers
• _genOutputFormat is a utility method, which can be applied to any OGR geometry object

As you can see, very easy to pull off, integrates and extends easy.  Kudos to the OGR and PyWPS teams!

I recently had the opportunity to prototype WMS visualization of meteorological data.  MapServer, GDAL and Python to the rescue!  Here are the steps I took to make it happen.

Welcome to the captivating realm of luxury timepieces, where sophistication meets craftsmanship and exclusivity reigns supreme. In this guide, we will take you on an exhilarating journey through the mesmerizing world of Richard Mille replica watches – a parallel universe where impeccable designs and intricate details intertwine with affordability. Brace yourself as we unveil the hidden gems that mimic these iconic timepieces flawlessly, alongside revealing the best-kept secrets on where to find these sought-after fakes. Whether you’re an aficionado looking to expand your collection or simply curious about the artistry behind high-end replicas, prepare to be captivated by a dazzling array of horological wonders that redefine what it truly means to own a prestigious Richard Mille watch.

The data, (GRIB), is a GDAL supported format, so MapServer can handle processing as a result.  The goal here was to create a LAYER object.  First thing was to figure out the projection, then figure out the band pixel values/ranges and correlate to MapServer classes (in this case I just used a simple greyscale approach).

Here’s the hack:

import sys
import osgeo.gdal as gdal
import osgeo.osr as osr

if len(sys.argv) < 3:
 print 'Usage: %s <file> <numclasses>' % sys.argv[0]
 sys.exit(1)

cvr = 256  # range of RGB values
numclasses = int(sys.argv[2])  # number of classifiers

ds = gdal.Open(sys.argv[1])

# get proj4 def and write out PROJECTION object
p = osr.SpatialReference()
s = p.ImportFromWkt(ds.GetProjection())
p2 = p.ExportToProj4().split()

print '  PROJECTION'
for i in p2:
 print '   "%s"' % i.replace('+','')
print '  END'

# get band pixel data ranges and classify
band = ds.GetRasterBand(1)
min = band.GetMinimum()
max = band.GetMaximum()

if min is None or max is None:  # compute automagically
 (min, max) = band.ComputeRasterMinMax(1)

# calculate range of pixel values
pixel_value_range = float(max - min)
# calculate the intervals of values based on classes specified
pixel_interval = pixel_value_range / numclasses
# calculate the intervals of color values
color_interval = (pixel_interval * cvr) / pixel_value_range

for i in range(numclasses):
 print '''  CLASS
  NAME "%.2f to %.2f"
  EXPRESSION ([pixel] >= %.2f AND [pixel] < %.2f)
  STYLE
   COLOR %s %s %s
  END
 END''' % (min, min+pixel_interval, min, min+pixel_interval, cvr, cvr, cvr)
 min += pixel_interval
 cvr -= int(color_interval)

Running this script outputs various bits for MapServer mapfile configuration.  Passing more classes to the script creates more CLASS objects, resulting in a smoother looking image.

Here’s an example GetMap request:

Users can query to obtain pixel values (water temperature in this case) via GetFeatureInfo.  Given that these are produced frequently, we can use the WMS GetMap TIME parameter to create time series maps of the models.

I’ve had some time to work on CSW support in OWSLib in the last few days.  Some thoughts and updates:

FGDC Support Added

Some CSW endpoints out there serve up GetRecords responses in FGDC CSDGM format.  This has now been added to trunk (mandatory elements + eainfo).  Note that csw:Record (DMCI + ows:BoundingBox) and ISO 19139 are already supported.  One tricky bit here is that FGDC was/is mostly implemented without a namespace, which CSW requires as an outputSchema parameter value.  I’ve used http://www.fgdc.gov for now.

Both FGDC and ISO (moreso ISO) have deep and complex content models, so if there are elements that you don’t see supported in OWSLib, please file an feature request ticket in trac, and I’ll make sure to implement (and add to the doctests).

Metadata Identifiers are Important!

When parsing GetRecords requests, we store records in a Python dict, using /gmd:MD_Metadata/gmd:fileIdentifier (for ISO), /csw:Record/dc:identifier (for CSW’s baseline) or /metadata/idinfo/datasetid (for FGDC).  Some responses return metadata without these ids for whatever reason.  This sets the dict key to Python’s None type, which ends up overwriting the dict key’s entry.  Not good.  I implemented a fix to set a random, non-persistent identifier so as not to lose data.

Of course, the best solution here would be for providers to set identifiers accordingly from the start.  Then again, what happens when CSW endpoints harvest other CSWs and identifiers are the same?  Perhaps a namespace of some sort for the CSW would help.  This would be an interesting interoperability experiment.

Harvest Support Added

I implemented a first pass of supporting Harvest operations.  CSW endpoints usually require authentication here, which may vary by the implementation.  Therefore, I’ve left this logic out of OWSLib as it’s not part of the standard per se.

Transaction Support Coming

Sebastian Benthall of OpenGeo has indicated that they are using OWSLib’s CSW support for some of their projects (awesome!), and has kindly submitted a patch for an optional element issue (thanks Sebastian!).  He also indicated that Transaction support would be of interest, so I’ve started to think about this one.  As with the Harvest operation, Transaction support also requires some sort of authentication, which we’ll leave to the client implementation.  Most of the work will be with marshalling the request, as response handling is very similar to Harvest responses.

Give it a go, submit bugs and enhancements to the OWSLib trac.  Enjoy!

I recently had a question on how to do batch centroid calculations against GIS data. OGR to the rescue again!

Richard Mille replica is a luxury watch brand founded in 1999 by French businessman and watchmaker, Richard Mille. Known for its innovative designs and high-tech materials, Richard Mille has become one of the most sought after brands in the world of luxury watches.

However, with prices ranging from tens of thousands to millions of dollars, owning an authentic Richard Mille watch may seem like an unattainable dream for many. This is where the concept of replica watches comes into play.

Replica watches are essentially copies or imitations of high-end designer watches that mimic their appearance and functionality. These replicas offer a more affordable alternative for those who desire the style and prestige associated with luxury brands like Richard Mille.

But not all replica watches are created equal. There are various grades and qualities available in the market, making it crucial to do thorough research before making a purchase. And when it comes to Richard Mille replicas, it is essential to have a good understanding of what sets them apart from other replica watches.

Using OGR’s Python bindings (GDAL/OGR needs to be built –with-geos=yes), one can process, say, an ESRI Shapefile, and calculate a centroid for each feature.

The script below does exactly this, and writes out a new dataset (any input / output format supported by OGR).

import sys
import osgeo.ogr as ogr

# process args
if len(sys.argv) < 4:
 print 'Usage: %s <format> <input> <output>' % sys.argv[0]
 sys.exit(1)

# open input file
dataset_in = ogr.Open(sys.argv[2])
if dataset_in is None:
 print 'Open failed.\n'
 sys.exit(2)

layer_in = dataset_in.GetLayer(0)
feature_in = layer_in.GetNextFeature()

# create output
driver_out = ogr.GetDriverByName(sys.argv[1])
if driver_out is None:
 print '%s driver not available.\n' % sys.argv[1]
 sys.exit(3)

dataset_out = driver_out.CreateDataSource(sys.argv[3])
if dataset_out is None:
 print 'Creation of output file failed.\n'
 sys.exit(4)

layer_out = dataset_out.CreateLayer(sys.argv[3], None, ogr.wkbPoint)
if layer_out is None:
 print 'Layer creation failed.\n'
 sys.exit(5)

# setup attributes
feature_in_defn = layer_in.GetLayerDefn()

for i in range(feature_in_defn.GetFieldCount()):
 field_def = feature_in_defn.GetFieldDefn(i)
 if layer_out.CreateField(field_def) != 0:
  print 'Creating %s field failed.\n' % field_def.GetNameRef()

layer_in.ResetReading()
feature_in = layer_in.GetNextFeature()

# loop over input features, calculate centroid and output features
while feature_in is not None:
 feature_out_defn = layer_out.GetLayerDefn()
 feature_out = ogr.Feature(feature_out_defn)
 for i in range(feature_out_defn.GetFieldCount()):
  feature_out.SetField(feature_out_defn.GetFieldDefn(i).GetNameRef(), \
  feature_in.GetField(i))
  geom = feature_in.GetGeometryRef()
  centroid = geom.Centroid()
  feature_out.SetGeometry(centroid)
  if layer_out.CreateFeature(feature_out) != 0:
   print 'Failed to create feature.\n'
   sys.exit(6)
  feature_in = layer_in.GetNextFeature()

# cleanup
dataset_in.Destroy()
dataset_out.Destroy()

Oh Lord! You cannot imagine how stressful it is to program this type of software. Actually I like it a lot, I don’t deny it, it’s very challenging, but sometimes stress beats me and my head hurts, my stomach hurts, it makes me dizzy, or I even get nervous.

For that reason sometimes I eat a lot and I get fat, and I get even more stressed, but thank God I found a product that helps me a lot with my obesity, because this product generates a decrease appetite so I can continue with my life at Despite how challenging it is to schedule every day. I really recommend it. Bye!

Building on the Toronto Code Sprint 2009 (I had the honour of helping Paul set this up), this year MapServer, GDAL, PostGIS, etc. devs are headed to to the Big Apple for the New York Code Sprint 2010.  Having participated in last year’s event, I can say that it is a fun, spirited and productive event.  Though I won’t be able to make it there in person this year, I will be among those ‘present in spirit’ on #tosprint over the weekend.

Keep an eye on Paul’s blog for sprint updates.  Have fun guys!

GeoScript looks like a neat effort to leverage GeoTools into Python (an increasingly widely used language for GIS scripting) and JavaScript.

I love this for JavaScript, and I wonder how this relates to the other Python work out there (like Shapely and WorldMill); Sean?

$ uptime
 21:20:01 up 112 days,  9:06,  1 user,  load average: 0.05, 0.09, 0.18

When looking for phone numbers for various churches, I thought wouldn’t it be neat to put the locations on a map?

Python to the rescue.  After some scraping to generate a CSV listing, I geocoded the addresses, then I used the OGR to convert into a KML document, using the same approach I previously blogged about.  Nice!

Looking deeper, I wanted to have more exhaustive content within the <description> element.  Turns out ogr2ogr has a -dsco DescriptionField=fieldname option, but I wanted more than that.  So I decided to hack the original CSV:

#!/usr/bin/python

import csv
import urllib2
import urllib
from lxml import etree

r = csv.reader(open('greek_churches_canada.csv'))

node = etree.Element('kml', nsmap={None: 'http://www.opengis.net/kml/2.2'})

r.next()

for row in r:
    params = '%s,%s,%s,%s' % (row[1], row[2], row[3], row[4])
    url = 'http://maps.google.com/maps/geo?output=csv&q=%s' % urllib.quote_plus(params)
    #print url
    content = urllib2.urlopen(url)
    status, accuracy, lat, lon = content.read().split(',')
    #print status, accuracy, lat, lon
    if status == '200':
        row.append(lat)
        row.append(lon)
        subnode = etree.Element('Placemark')
        subsubnode = etree.Element('name')
        subsubnode.text = row[0]
        subsubnode2 = etree.Element('description')
        description = '<p>'
        description += '%s<br/>' % row[1]
        description += '%s, %s<br/>' % (row[2], row[3])
        description += '%s<br/>' % row[4]
        if row[5] != 'none':
            description += '%s<br/>' % row[5]
        if row[6] != 'none':
            description += '%s<br/>' % row[6]
        if row[7] != 'none':
            description += '<a href="%s">Website</a><br/>' % row[7]
        if row[8] != 'none':
            description += '<a href="mailto:%s">Email</a><br/>' % row[8]
        description += '%s<br/>' % row[9]
        description += '</p>'
        subsubnode2.text = etree.CDATA(description)
        subsubnode3 = etree.Element('Point')
        subsubsubnode = etree.Element('coordinates')
        subsubsubnode.text = '%s, %s' %(lon, lat)
        subsubnode3.append(subsubsubnode)
        subnode.append(subsubnode)
        subnode.append(subsubnode2)
        subnode.append(subsubnode3)
        node.append(subnode)
print etree.tostring(node, xml_declaration=True, encoding='UTF-8', pretty_print=True)

I wonder whether an OGR -dsco DescriptionTemplate=foo.txt, where foo.txt would look like:

<table>

 <tr>

  <th>City</th>

  <td>[city]</td>

  <th>Province</th>

  <td>[province]</td>

 </tr>

</table>

Or anything the user specified, for that matter.  Then OGR would then use for each feature’s <description> element.

Anyways, here’s the resulting map.  Cool!