tommy's scratchpad

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!

Great post by Dave here on his experience and suggestions / ideas on finding a job.  Upbeat, positive and encouraging.  Congratulations and good post Dave!

There are several ways to apply for a job, one is through virtual means, but the other where there are more options is the physical means, once you get a job the best thing you can do with that money is login pragmatic play and shoot yourself to heaven with the profits you can make.

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.