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"
    lyrobj = None
    wfs = None
    cache = False # do not cache by default

    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")
        layerobj = self.mapobj.getLayerByName(typename)

        crs = self.__getLayerCrs(layer.crsOptions)

        if not crs:
            crs = Crs("epsg:4326")

        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)

        # do not take anything behind max extent
        #bbox[0] = layerobj.extent.minx if layerobj.extent.minx >= bbox[0] else bbox[0]
        #bbox[1] = layerobj.extent.miny if layerobj.extent.miny >= bbox[1] else bbox[1]
        #bbox[2] = layerobj.extent.maxx if layerobj.extent.maxx <= bbox[2] else bbox[2]
        #bbox[3] = layerobj.extent.maxy if layerobj.extent.maxy <= bbox[3] else bbox[3]

        bbox = map(lambda x: round(float(x),4), bbox)

        # 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: round(float(x),4), open(bboxfile).read().split(","))

                    # compare bounding boxes
                    if self.cache and \
                      (storedbbox[0] <= bbox[0] and\
                       storedbbox[1] <= bbox[1] and \
                       storedbbox[2] >= bbox[2] and \
                       storedbbox[3] >= bbox[3]) or \
                       self.request.upper() != "GETMAP":
                            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 = None
            if self.capabilities.version == "2.0.0":
                feature = self.capabilities.getfeature(
                        typename=typename, 
                        bbox=bbox,
                        filter=fes,
                        propertyname=None)
            else:
                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()

        # check
        if self.request.upper() in  ["GETMAP"]:
            check = ogr.Open(outfn)
            if not check:
                self.__clear(datadir,typename)
            else:
                check.Destroy()

        # set layer connection
        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.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)
            extent = layer.boundingBox
            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(layer)

            lyrobj.setProjection(crs.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
        """
        
        
        featureFile = open(os.path.join(self.cachedir,"feature.gml"),"w")
        featureFile.write(self.capabilities.getfeature( layer.id, count=1).read())
        featureFile.close()
        
        ds = ogr.Open(featureFile.name)
        ogrlayer = ds.GetLayer()
        geomType = ogrlayer.GetGeomType()
        if geomType == 0: # unknown
            # brutal force way
            f = ogrlayer.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
        """

        if len(crss) > 0:
            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))