Proxy4OWS/wfs/__init__.py

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

from OWS import OWS
import mapscript
import cgi
from lxml import objectify
import urllib
import urlparse
import logging
from osgeo import ogr,gdal,osr
from owslib.crs import Crs
import os,sys
import re
import shutil
import tempfile

class WFS(OWS):

    service = "WFS"
    wfsns = "http://www.opengis.net/wfs"
    layerDefFile = None
    lyrobj = None
    wfs = None

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

    def dispatch(self):
        """Dispatch given request
        """
        request = mapscript.OWSRequest()
        request.loadParams()

        typename = request.getValueByName("layers")
        if typename:
            layer = self.capabilities.contents[typename]

        self.request=request.getValueByName("REQUEST")

        # if no 'map' parameter in URL, create new mapfile
        if not request.getValueByName("map"):
            logging.debug("Creating new mapfile")
            self.mapobj = self.makeMap()
        else:
            # there is 'map' parameter in URL and the file exists, load it
            logging.debug("Using existing mapfile %s" % request.getValueByName("map"))
            if os.path.isfile(request.getValueByName("map")):
                self.mapobj = mapscript.mapObj(request.getValueByName("map"))
            # there is 'map' parameter in URL BUT the file does not exist:
            # create
            else:
                self.mapobj = self.makeMap(request.getValueByName("map"))

            if self.mapobj:
                self.mapfilename = request.getValueByName("map")

        # download data
        if self.request.upper() in  ["GETMAP", "GETFEATUREINFO"]:
            dataFile = self.getData(request,typename,layer)
            if not ogr.Open(dataFile):
                dataFile = self.getData(request,typename,layer)
            if dataFile:
                #layer.data = dataFile.name
                pass


        mapscript.msIO_installStdoutToBuffer()
        res = self.mapobj.OWSDispatch(request)

        if mapscript.MS_DONE == res:
            raise OWSExceptions.NoValidRequest("No valid OWS Request")
        elif mapscript.MS_FAILURE == res:
            pass
            #raise OWSExceptions.RequestFailed("Request failed")


        print "Content-type: %s\n" % (mapscript.msIO_stripStdoutBufferContentType())
        # HACK HACK HACK
        # Fix missing namespace encoding for ArcIMS and similar layers
        if self.request.upper() == "GETFEATUREINFO" and\
                typename.find(":") > -1:
            print self.__fixNamespace(typename)
        else:
            print mapscript.msIO_getStdoutBufferBytes()

    def getData(self,request,typename,layer):
        """Download data from WFS server and store them to local harddrive
        """

        fes = request.getValueByName("fes")
        bbox = request.getValueByName("bbox").split(",")
        featureid = request.getValueByName("featureid")
        featureversion = request.getValueByName("featureversion")
        datadir = os.path.join(self.cachedir,"cache")

        layer = self.capabilities.contents[typename]
        crs = self.__getLayerCrs(layer.crsOptions)
        version = request.getValueByName("version")

        # get propper bbox for the WFS request
        if version == "1.3.0":
            requestCrs = request.getValueByName("crs")
        else:
            requestCrs = request.getValueByName("srs")
        requestCrs = Crs(requestCrs)

        bbox = self.__adjustBBox(bbox,requestCrs, crs,version)

        # clear dir
        outfn = None 

        if os.path.isdir(datadir):
            if os.path.isfile(os.path.join(datadir,"%s.gml"%typename)):
                # find out bbox of the data
                bboxfile = os.path.join(datadir,"%s.bbox"%typename)
                filterfile = os.path.join(datadir,"%s.filter"%typename)
                if os.path.isfile(bboxfile) and not os.path.isfile(filterfile)\
                        and not fes:
                    # convert "string" to [floats] as [minx,miny,maxx,maxy]
                    storedbbox = map(lambda x: float(x), open(bboxfile).read().split(","))
                    if round(storedbbox[0],4) <= round(bbox[0],4) and\
                    round(storedbbox[1],4) <= round(bbox[1],4) and \
                    round(storedbbox[2],4) >= round(bbox[2],4) and \
                    round(storedbbox[3],4) >= round(bbox[3],4):
                        logging.info(
                                "Using cached file for type [%s] with bbox [%f,%f,%f,%f], not downloading new data"%\
                                        (typename, bbox[0],bbox[1],bbox[2],bbox[3]))
                        # setting output file name
                        outfn = os.path.join(datadir,"%s.gml"%typename)
                    else:
                        # remove pre-cached file with only little area
                        logging.info("Removing pre-cached files %s.*"% typename)
                        self.__clear(datadir,typename)
        else:
            os.mkdir(os.path.join(self.cachedir,"cache"))

        # create new cache file, if it does not exist yet
        # download the data from the server
        if outfn == None:
            # clear, just to be sure
            self.__clear(datadir,typename)
                
            # create cached bbox
            outbbox = open(os.path.join(datadir,"%s.bbox"%typename),"w")
            outbbox.write("%f,%f,%f,%f"%(bbox[0],bbox[1],bbox[2],bbox[3]))
            outbbox.close()

            # create cached filter
            if fes:
                outfes = open(os.path.join(datadir,"%s.filter"%typename),"w")
                outfes.write(fes)
                outfes.close()

            # create chace file
            outfn= os.path.join(datadir,"%s.gml"%typename)
            #cachefile = open(os.path.join(datadir,"%s.gmlorig"%typename),"w")
            cachefile = open(outfn,"w")

            # download feature from WFS
            logging.debug("Downloading data [%s] from bbox [%f,%f,%f,%f]"%\
                            (typename, bbox[0],bbox[1],bbox[2],bbox[3]))
            bbox.append(crs.getcode())
            feature = self.capabilities.getfeature(
                        typename=typename, 
                        bbox=bbox,
                        filter=fes,
                        srsname=crs.getcode(),
                        propertyname=None)

            # features are here
            # now take OGR and convert downloaded data into OGR-readable
            # format
            # this will ommit some GML-specific issues (like namespaces
            # etc.)
            # data will be stored into target file name

            cachefile.write(feature.read())
            cachefile.close()
            #ds = ogr.Open(cachefile.name)
            #drv = ogr.GetDriverByName(ds.GetDriver().name)
            #out = drv.CopyDataSource(ds,outfn)
            #out.Destroy()

        # set layer connection
        layerobj = self.mapobj.getLayerByName(typename)
        layerobj.connection = os.path.abspath(outfn)

        # data are downloaded
        # make sure, they have propper axis order - x,y
        if self.capabilities.version == "1.1.0" and crs.axisorder == "yx":
            pass
            logging.debug("Setting GML_INVERT_AXIS_ORDER_IF_LAT_LONG variable to YES")
            gdal.SetConfigOption("GML_INVERT_AXIS_ORDER_IF_LAT_LONG","YES")
            gdal.SetConfigOption("GML_CONSIDER_EPSG_AS_URN","YES")

            # >>> from osgeo import ogr
            # >>> from osgeo import gdal
            # >>> gdal.SetConfigOption("GML_INVERT_AXIS_ORDER_IF_LAT_LONG","YES")
            # >>> gdal.SetConfigOption("GML_CONSIDER_EPSG_AS_URN","YES")
            # >>> ind = ogr.Open("data.xml")
            # >>> outd = ogr.GetDriverByName("GML")
            # >>> outf = outd.CopyDataSource(ind,"out4.xml")
            # >>> outf.Destroy()

        logging.info("Connection of [%s] set to %s" % (typename,layerobj.connection))
        return outfn

    def __adjustBBox(self,bbox,src,target,version):
        """ adjust bounding coordinates and axis order
        """

        # swap axis, if needed
        if version == "1.3.0" and src.axisorder == "yx":
            bbox[0],bbox[1] = bbox[1],bbox[0]
            bbox[2],bbox[3] = bbox[3],bbox[2]

        # coordinate transformation
        projsrc = osr.SpatialReference()
        projsrc.ImportFromEPSG(src.code)

        projtarget = osr.SpatialReference()
        projtarget.ImportFromEPSG(target.code)

        trans = osr.CoordinateTransformation(projsrc,projtarget)
        bbox[0],bbox[1],z = trans.TransformPoint(float(bbox[0]),float(bbox[1]))
        bbox[2],bbox[3],z = trans.TransformPoint(float(bbox[2]),float(bbox[3]))

        return bbox

    def setFilter(self,mapobj,request):
        """ Set WFS filter encoding 
        """
        # get the layer
        layerobj = mapobj.getLayerByName(request.getValueByName("layers"))


        # get the filter
        logging.debug("FES received from HSLayers: %s" % fes)

        # cut off the opening and closing <Filter> tag
        # - this is needed for mapserver
        #root = etree.XML(fes)
        #msFilter = ""
        #for child in root:
        #    msFilter += etree.tostring(child)
        #logging.debug("Setting the filter %s" % msFilter)

        # set the filter
        layerobj.setMetaData("wfs_filter",fes)
        #layerobj.setMetaData("wfs_filter",msFilter)

        # save the mapfile - debugging
        mapobj.save("mapfile.fes")

    def makeMap(self,mapfilename=None):

        mapobj = self.getMapObj(mapfilename)

        # mapobjects exists, do not do any new layers
        if mapobj.numlayers:
            return mapobj

        self.layerDefFile = self.createLayerDefinitionFile("wfs",
                os.path.join( os.path.dirname(__file__), "templates",'wfs.xml'))

        ds = ogr.Open(self.layerDefFile)

        self.setMapName(mapobj)

        # load mapfile SYMBOLs
        symbolPath = os.path.join( os.path.dirname(__file__), "symbols.txt")
        symbolsLoaded = mapobj.setSymbolSet(symbolPath)
        if symbolsLoaded == mapscript.MS_SUCCESS:
            logging.debug("Symbols loaded from %s",symbolPath)
        else:
            logging.debug("Error loading symbols from %s",symbolPath)

        logging.debug(self.capabilities.contents)

        srss = []

        for name in self.capabilities.contents:
                
            mapobj.setMetaData("wms_srs",self.config.get("MapServer","srs"))
            mapobj.setMetaData("wfs_srs",self.config.get("MapServer","srs"))
            mapobj.setMetaData("wfs_connectiontimeout","90")
            layer = self.capabilities.contents[name]

            logging.debug("Creating layer %s" % name)

            srss = srss+filter(lambda y: not y in srss,layer.crsOptions)

            lyrobj = mapscript.layerObj(mapobj)
            #lyrobj.name = name.replace(":","_")
            lyrobj.name = name
            lyrobj.title = layer.title
            if layer.title:
                lyrobj.setMetaData("wms_title",layer.title)
                lyrobj.setMetaData("wfs_title",layer.title)
            #if layer.abstract:
            #    lyrobj.setMetaData("ows_abstract",  layer.abstract)
            logging.debug("WFS version %s",self.capabilities.version)
            lyrobj.setMetaData("gml_include_items","all")
            lyrobj.setConnectionType(mapscript.MS_OGR,'')
            lyrobj.data = re.sub(r".*:","",name)
            crs = self.__getLayerCrs(layer.crsOptions)
            if ds:
                ogrLayer = ds.GetLayerByName(name)
                extent = self.getLayerExtent(layer,crs)
                if extent:
                    lyrobj.setMetaData("wms_extent","%s %s %s %s" % \
                            (extent[0],extent[1],extent[2],extent[3]))
                    lyrobj.setMetaData("wfs_extent","%s %s %s %s" % \
                            (extent[0],extent[1],extent[2],extent[3]))
                lyrobj.type = self._getLayerType(ogrLayer)
            else:
                mapobj.removeLayer(mapobj.numlayers-1)
                logging.debug("No ogrDataSource found")
                continue

            lyrobj.setProjection(crs.getcode())
            #lyrobj.setProjection(layer.crsOptions[0].getcode())

            lyrobj.dump = mapscript.MS_TRUE 
            lyrobj.template = "foo"
            cls = mapscript.classObj(lyrobj)
            style = mapscript.styleObj(cls)
            style.outlinecolor=mapscript.colorObj(134,81,0)
            style.color=mapscript.colorObj(238,153,0)
            style.size=5
            style.width=5

            if lyrobj.type == mapscript.MS_LAYER_POINT:
                style.symbol = 1

        ## overwrite already set SRSs
        #if len(srss) > 0:
        #    logging.debug("Overwriting SRS option")
        #    mapobj.setMetaData("wms_srs"," ".join(srss))

        self.saveMapfile(mapobj,mapfilename)
        return mapobj
    
    def getGeomName(self,geomname):

        if geomname.find("LINE") > -1:
            return mapscript.MS_LAYER_LINE
        elif geomname.find("POLYGON") > -1:
            return mapscript.MS_LAYER_POLYGON
        else:
            return mapscript.MS_LAYER_POINT

    def setMapName(self,mapobj):
        mapobj.name = self.config.get("MapServer","name")

        if self.capabilities.identification.title:
            mapobj.setMetaData("wms_title",self.capabilities.identification.title)
            mapobj.setMetaData("wfs_title",self.capabilities.identification.title)
        if self.capabilities.identification.abstract:
            mapobj.setMetaData("wms_abstract",self.capabilities.identification.abstract)
            mapobj.setMetaData("wfs_abstract",self.capabilities.identification.abstract)

    def _getLayerType(self,layer):
        """Returns MS layer type based on ogr.Layer.GetGeomType

        with ogr:

            wkbGeometryCollection = 7
            wkbGeometryCollection25D = -2147483641
            wkbLineString = 2
            wkbLineString25D = -2147483646
            wkbLinearRing = 101
            wkbMultiLineString = 5
            wkbMultiLineString25D = -2147483643
            wkbMultiPoint = 4
            wkbMultiPoint25D = -2147483644
            wkbMultiPolygon = 6
            wkbMultiPolygon25D = -2147483642
            wkbNDR = 1
            wkbNone = 100
            wkbPoint = 1
            wkbPoint25D = -2147483647
            wkbPolygon = 3
            wkbPolygon25D = -2147483645
            wkbUnknown = 0
        """

        geomType = layer.GetGeomType()
        if geomType == 0: # unknown
            # brutal force way
            f = layer.GetNextFeature()
            if f:
                gr = f.GetGeometryRef()
                geomType = gr.GetGeometryType()

        if geomType in [ogr.wkbPolygon,
                        ogr.wkbMultiPolygon,
                        ogr.wkbLinearRing]:
               return mapscript.MS_LAYER_POLYGON
        elif geomType in [ogr.wkbLineString,
                          ogr.wkbMultiLineString]:
               return mapscript.MS_LAYER_LINE
        else:
               return mapscript.MS_LAYER_POINT

    def __getLayerCrs(self,crss):
        """
        Returns bests (non-degree) coordinate system of the layer, which is
        available.

        Ofcourse, sometimes, there is no other option, there EPSG:4326, but
        take somethign else, if you can
        """
        for crs in crss:
            if crs.getcode() == "EPSG:4326":
                return crs
        return crss[0]

    def __clear(self,datadir,typename):
        """Remove all cached files with following typename
        """

        for i in os.listdir(datadir):
            if i.find(typename) == 0:
                logging.debug("Removing pre-cached file %s"% i)
                os.remove(os.path.join(datadir,i))

    def __fixNamespace(self,typename):
        """Fix missing namespace in output XML
        """
        content = mapscript.msIO_getStdoutBufferBytes()
        
        # look after something like
        # <namespace:layerName_layer> and add namespace definition
        
        if typename.find(":") == -1:
            return content

        (prefix,name) = typename.split(":")
        namespace = self.capabilities._capabilities.nsmap[prefix]
        return content.replace("<%s:%s%s>"%(prefix,name,"_layer"),
                               "<%s:%s%s xmlns:%s=\"%s\">"%(prefix,name,"_layer",prefix,namespace))