Proxy4OWS/wcs/__init__.py

#!/usr/bin/env python
# coding=utf-8

from OWS import OWS
import mapscript
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
import os
import logging
import re

class WCS(OWS):

    service = "WCS"
    request = '' # GetCapabilies, ...
    parser = None
    version = None

    def __init__(self,url=None,qstring=None,configFiles=None):
        OWS.__init__(self,url,qstring)

    def makeMap(self,mapfilename=None):
        """Create mapscript.mapObj"""

        lyrobj = None

        self.version = self.capabilities.version

        mapobj = self.getMapObj(mapfilename)

        #
        # for each layer from the WCS Capabilities file, transform it to
        # mapscript.layerObj
        #
        for name in self.capabilities.contents:
            layer = self.capabilities.contents[name]
            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
                if layer.title:
                    lyrobj.title = layer.title
                    lyrobj.setMetaData("wms_title",layer.title)
                #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"))
                extent = None
                if layer.crsOptions:
                    lyrobj.setProjection(layer.crsOptions[0])
                    extent = self.getLayerExtent(layer,layer.crsOptions[0])
                else:
                    sr = osr.SpatialReference()
                    sr.ImportFromWkt(ds.GetProjection())
                    if sr.AutoIdentifyEPSG() == 0:
                        epsg = "{code}:{value}".format(code=sr.GetAuthorityName("PROJCS"),value=sr.GetAuthorityCode("PROJCS"))
                        extent = self.getLayerExtent(layer,epsg)
                        lyrobj.setProjection(epsg)
                    else:
                        extent = self.getLayerExtent(layer,wkt=ds.GetProjection())
                        lyrobj.setProjection(sr.ExportToProj4())
                lyrobj.setMetaData("wms_extent","%s %s %s %s" % \
                        (extent[0],extent[1],extent[2],extent[3]))
                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("Layer %s: %s, skipping." %\
                                (name,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")