Geografic-Information-System/gui/wxpython/core/ws.py

"""
@package core.ws

@brief Fetching and preparation of web service data for rendering.

Note: Currently only WMS is implemented.

Classes:
 - ws::RenderWMSMgr
 - ws::GDALRasterMerger

(C) 2012 by the GRASS Development Team

This program is free software under the GNU General Public License
(>=v2). Read the file COPYING that comes with GRASS for details.

@author Stepan Turek <stepan.turek seznam.cz> (mentor: Martin Landa)
"""
import sys
import copy

import wx
from wx.lib.newevent import NewEvent

from grass.script.utils import try_remove
from grass.script import core as grass

from core          import utils
from core.debug    import Debug

from core.gconsole import CmdThread, GStderr, EVT_CMD_DONE, EVT_CMD_OUTPUT
from core.utils import _

try:
    haveGdal = True
    from osgeo import gdal
    from osgeo import gdalconst 
except ImportError:
    haveGdal = False

from grass.pydispatch.signal import Signal


class RenderWMSMgr(wx.EvtHandler):
    """Fetch and prepare WMS data for rendering.
    """
    def __init__(self, layer, mapfile, maskfile):
        if not haveGdal:
            sys.stderr.write(_("Unable to load GDAL Python bindings.\n"\
                               "WMS layers can not be displayed without the bindings.\n"))

        self.layer = layer

        wx.EvtHandler.__init__(self)

        # thread for d.wms commands
        self.thread = CmdThread(self)
        self.Bind(EVT_CMD_DONE, self.OnDataFetched)

        self.downloading = False
        self.renderedRegion = None
        self.updateMap = True
        self.fetched_data_cmd = None

        self.cmdStdErr = GStderr(self)

        self.mapfile = mapfile
        self.maskfile = maskfile
        self.tempMap = grass.tempfile()
        self.dstSize = {}
 
        self.Bind(EVT_CMD_OUTPUT, self.OnCmdOutput)
        
        self.dataFetched = Signal('RenderWMSMgr.dataFetched')
        self.updateProgress = Signal('RenderWMSMgr.updateProgress')

    def __del__(self):
        try_remove(self.tempMap)

    def Render(self, cmd, env):
        """If it is needed, download missing WMS data.

        .. todo::
            lmgr deletes mapfile and maskfile when order of layers
            was changed (drag and drop) - if deleted, fetch data again
        """
        if not haveGdal:
            return

        env = copy.copy(env)
        self.dstSize['cols'] = int(env["GRASS_WIDTH"])
        self.dstSize['rows'] = int(env["GRASS_HEIGHT"])

        region = self._getRegionDict(env)
        self._fitAspect(region, self.dstSize)

        self.updateMap = True
        fetchData = False
        zoomChanged = False

        if self.renderedRegion is None or \
           cmd != self.fetched_data_cmd:
            fetchData = True
        else:
            for c in ['north', 'south', 'east', 'west']:
                if self.renderedRegion and \
                   region[c] != self.renderedRegion[c]:
                    fetchData = True
                    break

            for c in ['e-w resol', 'n-s resol']:
                if self.renderedRegion and \
                    region[c] != self.renderedRegion[c]:
                    zoomChanged = True
                    break

        if fetchData:
            self.fetched_data_cmd = None
            self.renderedRegion = region

            try_remove(self.mapfile)
            try_remove(self.tempMap)

            self.currentPid = self.thread.GetId()
            self.thread.abort()
            self.downloading = True

            self.fetching_cmd = cmd
            cmdList = utils.CmdTupleToList(cmd)

            if Debug.GetLevel() < 3:
                cmdList.append('--quiet')

            env["GRASS_PNGFILE"] = self.tempMap
            env["GRASS_REGION"] = self._createRegionStr(region)

            self.thread.RunCmd(cmdList, env=env, stderr=self.cmdStdErr)

    def OnCmdOutput(self, event):
        """Print cmd output according to debug level.
        """
        if Debug.GetLevel() == 0:
            if event.type == 'error':
                sys.stderr.write(event.text)
                sys.stderr.flush()
        else:
            Debug.msg(1, event.text)

    def OnDataFetched(self, event):
        """Fetch data
        """
        if event.pid != self.currentPid:
            return
        self.downloading = False
        if not self.updateMap:
            self.updateProgress.emit(layer=self.layer)
            self.renderedRegion = None
            self.fetched_data_cmd = None
            return

        self.mapMerger = GDALRasterMerger(targetFile = self.mapfile, region = self.renderedRegion,
                                          bandsNum = 3, gdalDriver = 'PNM', fillValue = 0)
        self.mapMerger.AddRasterBands(self.tempMap, {1 : 1, 2 : 2, 3 : 3})
        del self.mapMerger

        self.maskMerger = GDALRasterMerger(targetFile = self.maskfile, region = self.renderedRegion,
                                           bandsNum = 1, gdalDriver = 'PNM', fillValue = 0)
        #{4 : 1} alpha channel (4) to first and only channel (1) in mask
        self.maskMerger.AddRasterBands(self.tempMap, {4 : 1}) 
        del self.maskMerger

        self.fetched_data_cmd = self.fetching_cmd

        self.dataFetched.emit()

    def _getRegionDict(self, env):
        """Parse string from GRASS_REGION env variable into dict.
        """
        region = {}
        parsedRegion = env["GRASS_REGION"].split(';')
        for r in parsedRegion:
            r = r.split(':')
            r[0] = r[0].strip()
            if len(r) < 2:
                continue
            try:
                if r[0] in ['cols', 'rows']:
                    region[r[0]] = int(r[1])
                else:
                    region[r[0]] = float(r[1])
            except ValueError:
                region[r[0]] = r[1]

        return region

    def _createRegionStr(self, region):
        """Create string for GRASS_REGION env variable from  dict created by _getRegionDict.
        """
        regionStr = ''
        for k, v in region.iteritems():
            item = k + ': ' + str(v)
            if regionStr:
                regionStr += '; '
            regionStr += item

        return regionStr

    def IsDownloading(self):
        """Is it downloading any data from WMS server? 
        """
        return self.downloading

    def _fitAspect(self, region, size):
        """Compute region parameters to have direction independent resolution.
        """
        if region['n-s resol'] > region['e-w resol']:
            delta = region['n-s resol'] * size['cols'] / 2

            center = (region['east'] - region['west'])/2

            region['east'] = center + delta + region['west']
            region['west'] = center - delta + region['west']
            region['e-w resol'] = region['n-s resol']

        else:
            delta = region['e-w resol'] * size['rows'] / 2

            center = (region['north'] - region['south'])/2 

            region['north'] = center + delta + region['south']
            region['south'] = center - delta + region['south']
            region['n-s resol'] = region['e-w resol']

    def Abort(self):
        """Abort process"""
        self.updateMap = False
        self.thread.abort(abortall = True)        


class GDALRasterMerger:
    """Merge rasters.

        Based on gdal_merge.py utility.
    """
    def __init__(self, targetFile, region, bandsNum, gdalDriver, fillValue = None):
        """Create raster for merging.
        """
        self.gdalDrvType = gdalDriver

        nsRes = (region['south'] - region['north']) / region['rows']
        ewRes = (region['east'] - region['west']) / region['cols']

        self.tGeotransform = [region['west'], ewRes, 0, region['north'], 0, nsRes]

        self.tUlx, self.tUly, self.tLrx, self.tLry = self._getCorners(self.tGeotransform, region)

        driver = gdal.GetDriverByName(self.gdalDrvType)
        self.tDataset = driver.Create(targetFile, region['cols'], region['rows'], bandsNum,  gdal.GDT_Byte)

        if fillValue is not None:
            # fill raster bands with a constant value
            for iBand in range(1, self.tDataset.RasterCount + 1):
                self.tDataset.GetRasterBand(iBand).Fill(fillValue)

    def AddRasterBands(self, sourceFile, sTBands):
        """Add raster bands from sourceFile into the merging raster.
        """
        sDataset = gdal.Open(sourceFile, gdal.GA_ReadOnly) 
        if sDataset is None:
            return

        sGeotransform = sDataset.GetGeoTransform()

        sSize = {
                    'rows' :  sDataset.RasterYSize,
                    'cols' :  sDataset.RasterXSize
                 }

        sUlx, sUly, sLrx, sLry = self._getCorners(sGeotransform, sSize)

        # figure out intersection region
        tIntsctUlx = max(self.tUlx,sUlx)
        tIntsctLrx = min(self.tLrx,sLrx)
        if self.tGeotransform[5] < 0:
            tIntsctUly = min(self.tUly,sUly)
            tIntsctLry = max(self.tLry,sLry)
        else:
            tIntsctUly = max(self.tUly,sUly)
            tIntsctLry = min(self.tLry,sLry)
        
        # do they even intersect?
        if tIntsctUlx >= tIntsctLrx:
            return
        if self.tGeotransform[5] < 0 and tIntsctUly <= tIntsctLry:
            return
        if self.tGeotransform[5] > 0 and tIntsctUly >= tIntsctLry:
            return


        # compute target window in pixel coordinates.
        tXoff = int((tIntsctUlx - self.tGeotransform[0]) / self.tGeotransform[1] + 0.1)
        tYoff = int((tIntsctUly - self.tGeotransform[3]) / self.tGeotransform[5] + 0.1)
        tXsize = int((tIntsctLrx - self.tGeotransform[0])/self.tGeotransform[1] + 0.5) - tXoff
        tYsize = int((tIntsctLry - self.tGeotransform[3])/self.tGeotransform[5] + 0.5) - tYoff

        if tXsize < 1 or tYsize < 1:
            return

        # Compute source window in pixel coordinates.
        sXoff = int((tIntsctUlx - sGeotransform[0]) / sGeotransform[1])
        sYoff = int((tIntsctUly - sGeotransform[3]) / sGeotransform[5])
        sXsize = int((tIntsctLrx - sGeotransform[0]) / sGeotransform[1] + 0.5) - sXoff
        sYsize = int((tIntsctLry - sGeotransform[3]) / sGeotransform[5] + 0.5) - sYoff

        if sXsize < 1 or sYsize < 1:
            return

        for sBandNnum, tBandNum in sTBands.iteritems():
            bandData = sDataset.GetRasterBand(sBandNnum).ReadRaster(sXoff, sYoff, sXsize,
                                                                    sYsize, tXsize, tYsize, gdal.GDT_Byte)
            self.tDataset.GetRasterBand(tBandNum).WriteRaster(tXoff, tYoff, tXsize, tYsize, bandData, 
                                                              tXsize, tYsize, gdal.GDT_Byte)

    def _getCorners(self, geoTrans, size):

        ulx = geoTrans[0]
        uly = geoTrans[3]
        lrx = geoTrans[0] + size['cols'] * geoTrans[1]
        lry = geoTrans[3] + size['rows'] * geoTrans[5]

        return ulx, uly, lrx, lry 

    def SetGeorefAndProj(self):
        """Set georeference and projection to target file
        """
        projection = grass.read_command('g.proj', 
                                        flags = 'wf')
        self.tDataset.SetProjection(projection)

        self.tDataset.SetGeoTransform(self.tGeotransform)

    def __del__(self):
        self.tDataset = None