added time to wcs

OWS.py

@@ -196,6 +196,7 @@ class OWS:
         mapobj.setProjection("init=epsg:4326")
         mapobj.setSize(500,500)
         mapobj.setExtent(-180,-90,90,180)
+        mapobj.shapepath = self.cachedir
 
         logging.debug("Setting ERRORFILE to %s"%self.config.get("MapServer","errorfile"))
         if not os.path.exists(self.config.get("MapServer","errorfile")):
@@ -241,13 +242,23 @@ class OWS:
 
         return layerurl
 
-    def createLayerDefinitionFile(self, name, templatefile):
+    def createLayerDefinitionFile(self, name,
+            templatefile,time=None,target=None):
 
             layerurl = self.getLayerUrl()
-            defFileName = os.path.join(self.cachedir,'%s.%s'%(name,self.service.lower()))
+            defFileName = None
+            if target:
+                defFileName = target
+            else:
+                defFileName = os.path.join(self.cachedir,'%s.%s'%(name,self.service.lower()))
+            if time:
+                time = "<DefaultTime>"+time+"</DefaultTime>"
+            else:
+                time = ""
+
             open(defFileName,'w').write(
                     Template(open(templatefile).read()).substitute( 
-                            dict(url= layerurl, name=name)))
+                            dict(url= layerurl, name=name,time=time)))
 
             logging.debug("Created %s layer definition file" % defFileName)
 

config.cfg-template

@@ -5,7 +5,7 @@ logging=DEBUG
 [MapServer]
 name=owsproxy
 tempdir=/tmp/
-errorfile=/tmp/mapserv.log
+errorfile=stderr
 imagepath=/tmp/mapserv/
 onlineresource=http://localhost/cgi-bin/owsproxy.cgi
 srs=EPSG:4326 EPSG:102067 EPSG:900913 EPSG:3035

wcs/__init__.py

@@ -5,9 +5,11 @@ from OWS import OWS
 import mapscript
 from osgeo import gdal
 from osgeo import osr
+from osgeo import ogr
 from string import Template
 import os
 import logging
+import re
 
 class WCS(OWS):
 
@@ -21,25 +23,33 @@ class WCS(OWS):
         OWS.__init__(self,url,qstring)
 
     def makeMap(self,mapfilename=None):
+        """Create mapscript.mapObj"""
 
         lyrobj = None
 
         self.version = self.capabilities.attrib["version"]
         self.wcsns = self.capabilities.nsmap[self.capabilities.prefix]
 
-
         mapobj = self.getMapObj(mapfilename)
 
+        #
+        # for each layer from the WCS Capabilities file, transform it to
+        # mapscript.layerObj
+        #
+        logging.info(self.capabilities["Contents"].getchildren())
         for layer in self.capabilities["{%s}Contents" % (self.wcsns)].getchildren():
             name =  layer.Identifier.text
             logging.debug("Creating layer %s" % name)
 
-
+            # create layer definition file
+            # see http://www.gdal.org/frmt_wcs.html
             layerDefFile = self.createLayerDefinitionFile(name,
                     os.path.join( os.path.dirname(__file__), "templates",'wcs.xml'))
 
+            # create the layer
             ds = gdal.Open(layerDefFile)
             try:
+
                 lyrobj = mapscript.layerObj(mapobj)
                 lyrobj.name = name
                 lyrobj.data = layerDefFile
@@ -47,18 +57,81 @@ class WCS(OWS):
                 lyrobj.setMetaData("wms_title",layer["{%s}%s" % (self.owsns11,"Title")].text)
                 if "{%s}%s" % (self.owsns11,"Abstract") in layer:
                     lyrobj.setMetaData("ows_abstract",  layer["{%s}%s" % (self.owsns11,"Abstract")].text)
+                lyrobj.setMetaData("wms_srs",self.config.get("MapServer","srs"))
                 sr = osr.SpatialReference()
                 sr.ImportFromWkt(ds.GetProjection())
                 lyrobj.setProjection(sr.ExportToProj4())
                 lyrobj.type = mapscript.MS_LAYER_RASTER
                 lyrobj.dump = mapscript.MS_TRUE
                 lyrobj.template = "foo"
+                self.getTime(ds,lyrobj)
                 cls = mapscript.classObj(lyrobj)
                 mapscript.styleObj(cls)
             except Exception,e:
                 mapobj.removeLayer(lyrobj.index)
-                logging.warning("Could not GetProjection of %s: %s, skipping." %\
-                                (layer.Identifier.text,e))
+                logging.warning("Layer %s: %s, skipping." %\
+                                (layer.Identifier.text,e.message))
 
         self.saveMapfile(mapobj,mapfilename)
         return mapobj
+
+    def getTime(self,ds,lyrobj):
+        """Generate time index file, set approripate layer metadata"""
+
+        # make time index file
+        drv = ogr.GetDriverByName("ESRI Shapefile")
+        tindex = '%s.%s.%s'%(lyrobj.name,self.service.lower(),"tindex")
+        logging.info("Creating time index definition file %s.shp in %s" % (tindex,self.cachedir))
+        tds = drv.CreateDataSource(os.path.join(self.cachedir,tindex))
+        lyr = tds.CreateLayer(tindex)
+        tile_field = ogr.FieldDefn("location", ogr.OFTString )
+        tile_field.SetWidth( 256 )
+        lyr.CreateField(tile_field)
+
+        time_field = ogr.FieldDefn("time", ogr.OFTString )
+        time_field.SetWidth( 256 )
+        lyr.CreateField(time_field)
+
+        timeMetadata = ""
+        for sds in ds.GetSubDatasets():
+            if sds[0].find('time') > -1:
+                # crate mapfile metadata
+                sdsname = sds[0].split(",")[0]
+                sdsvalue = re.findall(r'".*"',sdsname)[0].replace('"',"")
+                timeMetadata = timeMetadata+","+sdsvalue
+
+                # create OGR tindex vector feature
+                feature = ogr.Feature(lyr.GetLayerDefn())
+
+                # feature geometry
+                geotransform = ds.GetGeoTransform()
+                w = ds.RasterXSize
+                h = ds.RasterYSize
+                bbox = {"minx":geotransform[0],"miny":geotransform[5]*h+geotransform[3],"maxx":geotransform[1]*w+geotransform[0],"maxy":geotransform[3]}
+                poly = ogr.CreateGeometryFromWkt("""POLYGON((%(minx)s %(miny)s, %(maxx)s %(miny)s, %(maxx)s %(maxy)s, %(minx)s %(maxy)s, %(minx)s %(miny)s))""" % bbox)
+                feature.SetGeometry(poly)
+
+                # set time attribute
+                feature.SetField("time","%s" % (sdsvalue))
+                feature.SetField("location",sds[0])
+                lyr.CreateFeature(feature)
+
+        if timeMetadata.endswith(","):
+            timeMetadata = re.sub(r",$","",timeMetadata)
+        if timeMetadata.startswith(","):
+            timeMetadata = re.sub(r"^,","",timeMetadata)
+
+        if timeMetadata:
+            logging.info("Setting %s to %s layer" % (timeMetadata, lyrobj.name))
+
+            lyrobj.setMetaData("wms_timeextent", timeMetadata) 
+            lyrobj.setMetaData("wms_timeitem", "time") 
+
+            # NOTE: lyrobj.data will not be the wcs definition file
+            # anymore, BUT the time index shapefile
+            logging.info("Setting layer data to %s" %
+                    (os.path.join(tindex,tindex+".shp")))
+            lyrobj.data = None
+            lyrobj.tileindex = os.path.join(tindex,tindex+".shp")
+        else:
+            logging.info("No time subset found")

wcs/templates/wcs.xml

@@ -2,4 +2,5 @@
   <ServiceURL>$url</ServiceURL>
   <CoverageName>$name</CoverageName>
   <Timeout>5000</Timeout>
+  $time
 </WCS_GDAL>