Archive for web

Displaying GRIB data with MapServer

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.

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]

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)
   COLOR %s %s %s
 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:

Meteorological data in GRIB format via MapServer WMS

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.

OWSLib CSW Updates and Implementation Thoughts

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 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!

easy CSW with eXcat

If you have existing metadata and simply want a CSW interface as a means to search and discover your geospatial metadata, eXcat provides a simple solution.  Following the installation steps, it’s quite simple to populate your CSW:

$ cd excat/csw/WEB-INF/harvest
$ tar zxf my_metadata_files.tgz
$ for i in *.xml
> do
> lwp-download "http://localhost/excat/csw?request=Harvest&service=CSW&\
> version=2.0.2\&namespace=xmlns(csw=\
> source=$i&resourceFormat=application/xml\
> &resourceType="
> done

That’s pretty much it.  Lightweight simple approach, particularly for those who already have metadata management tools in place and need CSW.

/me thinks it would be nice to have a Python port of this for those in Java-less environments (why does it seem most CSW server implementations are in Java?)


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?

Python, KML and Parishes

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:


import csv
import urllib2
import urllib
from lxml import etree

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

node = etree.Element('kml', nsmap={None: ''})

for row in r:
    params = '%s,%s,%s,%s' % (row[1], row[2], row[3], row[4])
    url = '' % urllib.quote_plus(params)
    #print url
    content = urllib2.urlopen(url)
    status, accuracy, lat, lon =',')
    #print status, accuracy, lat, lon
    if status == '200':
        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)
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:


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!

Friday Metadata Thoughts

Not the most exciting topic, but I’ve found myself knee deep in metadata standards as they pertain to CSW in the last couple of weeks.

I’ve made some recommendations in the past for OWS metadata, which have helped in established publishing requirements for cataloguing.

Starting to look at ISO metadata (data, service) makes you quickly realize the pros and cons which come with making a standard flexible and exhaustive.  Let’s take 19139; almost everything in the schema is optional.  I think this is where profiles (such as ISO North American Profile) start to become especially important.

19119 is in the same boat.  Aside: then you start to wonder about the overlap between 19119 and OWS Capabilities metadata.  Wouldn’t it be nice if GetCapabiilties returned a 19119 document instead?  Which could plop nicely in a CSW query response as well. Oh wait, it already does.  But then try to validate the document instance.  You’ll find that OGC CSW and ISO use different versions of GML (follow the refs in the .xsd’s you’ll see them soon enough), yet apply them to the same namespace.  So validation fails.  Harmonization required!

Having said this, this is very complicated metadata which can be addressed by intelligent tools.  Tools that:

  • integrate with GIS systems which can automagically populate by:
    • fetching spatial extents
    • fetching reference system definitions
    • establish hierarchy (this would be tough as it would be tied to the data management of the system)
    • fetch contact information from a given user profile (how about getting this from the network’s email / global address book against the logged in user?)

Then again, what about keeping it simple and mainstream friendly?  The toughest part is metadata creation, so let’s make it as easy as possible to do so!

How do your activities try to make metadata easier to create?

new stuff in OWSLib

I’ve been spending alot of time lately doing a CSW client library in python, which was committed today to OWSLib.  CSW requests can be tricky to construct correctly, so this contribution attempts to provide an easy enough entry point to querying OGC Catalogues.

At this point, you can query your favourite CSW server with:

>>> from owslib import csw
>>> c = csw.request('')
>>> c.GetCapabilities() # constructs XML request, in c.request
>>> c.fetch() # HTTP call.  Result in c.response
>>> c.GetRecords('dataset','birds',[-152,42,-52,84]) # birds datasets in Canada       
>>> c.fetch()
>>> c.GetRecords('service','frog') # look for services with frogs, anywhere
>>> c.fetch()

That’s pretty much all there is to it.  There’s also support for DescribeRecord, GetRecordById and GetDomain for the adventurous.

I hope this will be a valuable addition.  Because CSW uses Filter, I broke things out into a module per standard, so that other code can reuse, say, filter for request building and response parsing.  A colleague is using this functionality to write a QGIS CSW search plugin.

My next goal will be to put in some response handling.  This will be tricky given the various outputSchema’s a given CSW advertises.  For now, I will concentrate on the default csw:Record (a glorified Dublin Core with ows:BoundingBox).

So try it out;  comments, feedback and suggestions would be most valued.

Oh ya, thank you etree!

MapServer 5.4.0 released

Announced yesterday, this release closes 92 bugs, and adds some new goodies.

Next stop: MapServer 6.0

Creating sitemap files for GeoNetwork

Sitemaps are a valuable way to index your content for web crawlers.  GeoNetwork is a great tool for metadata management and a portal environment for discovery.  I wanted to push out all metadata resources out as a sitemap so that content can be found by web crawlers.  Python to the rescue:

import MySQLdb
# connect to db
db=MySQLdb.connection(host='', user='foo',passwd='foo',db='geonetwork')
# print out XML header
print """<?xml version="1.0" encoding="UTF-8"?>

# fetch all metadata
db.query("""select id, schemaId, changeDate from Metadata where isTemplate = 'n'""")
r = db.store_result()

for row in r.fetch_row(0): # write out a url element
    if row[1] == 'fgdc-std':
        url = ''
    if row[1] == 'iso19139':
        url = ''
    print """ <url>
 </url>""" % (url, row[0], row[2], row[1])
print '</urlset>'

Done!  It would be great if this were an out-of-the-box feature of GeoNetwork.

MapServer Disaster: you have got to be kidding me

I’m beyond words at this point.

Modified: 1 April 2009 15:27:24 EST