""" @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 (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