[GRASS-SVN] r51352 - in grass-addons/grass7/raster: . r.in.wms2

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Apr 10 08:21:59 EDT 2012


Author: turek
Date: 2012-04-10 05:21:59 -0700 (Tue, 10 Apr 2012)
New Revision: 51352

Added:
   grass-addons/grass7/raster/r.in.wms2/
   grass-addons/grass7/raster/r.in.wms2/Makefile
   grass-addons/grass7/raster/r.in.wms2/r.in.wms2.py
   grass-addons/grass7/raster/r.in.wms2/wms_base.py
   grass-addons/grass7/raster/r.in.wms2/wms_drv.py
   grass-addons/grass7/raster/r.in.wms2/wms_gdal_drv.py
Log:
new module r.in.wms2 to download data from WMS server


Added: grass-addons/grass7/raster/r.in.wms2/Makefile
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/Makefile	                        (rev 0)
+++ grass-addons/grass7/raster/r.in.wms2/Makefile	2012-04-10 12:21:59 UTC (rev 51352)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.in.wms2
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

Added: grass-addons/grass7/raster/r.in.wms2/r.in.wms2.py
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/r.in.wms2.py	                        (rev 0)
+++ grass-addons/grass7/raster/r.in.wms2/r.in.wms2.py	2012-04-10 12:21:59 UTC (rev 51352)
@@ -0,0 +1,161 @@
+#!/usr/bin/env python
+"""
+MODULE:    r.in.wms2
+
+AUTHOR(S): Stepan Turek <stepan.turek AT seznam.cz>
+
+PURPOSE:   Downloads and imports data from WMS server.
+
+COPYRIGHT: (C) 2012 Stepan Turek, and 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.
+"""
+
+#%module
+#% description: Downloads and imports data from WMS servers.
+#% keywords: raster
+#% keywords: import
+#% keywords: wms
+#%end
+
+#%option
+#% key: mapserver
+#% type: string
+#% description:URL of WMS server 
+#% required: yes
+#% answer: http://wms.cuzk.cz/wms.asp
+#%end
+
+#%option
+#% key: layers
+#% type: string
+#% description: Layers to request from map server
+#% multiple: yes
+#% required: yes
+#% answer: prehledka_kraju-linie
+#%end
+
+#%option G_OPT_R_OUTPUT
+#% description: Name for output raster map
+#%end
+
+#%option
+#% key: srs
+#% type: integer
+#% description: EPSG number of source projection for request 
+#% answer:4326
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: region
+#% type: string
+#% description: Named region to request data for. Current region used if omitted
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: wms_version
+#% type:string
+#% description:WMS standard
+#% options:1.1.1,1.3.0
+#% answer:1.1.1
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: format
+#% type: string
+#% description: Image format requested from the server
+#% options: geotiff,tiff,jpeg,gif,png
+#% answer: geotiff
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: method
+#% type: string
+#% description: Reprojection method to use
+#% options:near,bilinear,cubic,cubicspline
+#% answer:near
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: maxcols
+#% type:integer
+#% description: Maximum columns to request at a time
+#% answer:400
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: maxrows
+#% type: integer
+#% description: Maximum rows to request at a time
+#% answer: 300
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: urlparams
+#% type:string
+#% description: Addition query parameters for server (only with 'd' flag)
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: styles
+#% type: string
+#% description: Styles to request from map server
+#% multiple: yes
+#% guisection: Map style
+#%end
+
+#%option
+#% key: bgcolor
+#% type: string
+#% description: Color of map background (only with 'd' flag)
+#% guisection: Map style
+#%end
+
+#%flag
+#% key: t
+#% description: Transparent background
+#% guisection: Map style
+#%end
+
+#%flag
+#% key: c
+#% description: Get capabilities
+#% guisection: Request properties
+#% suppress_required: yes
+#%end
+
+#%flag
+#% key: d
+#% description: Do not use GDAL WMS driver
+#%end
+
+import sys
+
+import grass.script as grass
+
+def main():
+    if flags['d']:
+        grass.debug("Using own driver")
+        from wms_drv import WMSDrv
+        wms = WMSDrv(options, flags)
+    else:
+        grass.debug("Using GDAL WMS driver")
+        from wms_gdal_drv import WMSGdalDrv
+        wms = WMSGdalDrv(options, flags)
+
+
+    
+    return 0
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    sys.exit(main())


Property changes on: grass-addons/grass7/raster/r.in.wms2/r.in.wms2.py
___________________________________________________________________
Added: svn:executable
   + *
Added: svn:mime-type
   + text/x-python
Added: svn:eol-style
   + native

Added: grass-addons/grass7/raster/r.in.wms2/wms_base.py
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/wms_base.py	                        (rev 0)
+++ grass-addons/grass7/raster/r.in.wms2/wms_base.py	2012-04-10 12:21:59 UTC (rev 51352)
@@ -0,0 +1,367 @@
+import os
+from   math import ceil
+
+import xml.etree.ElementTree as etree
+from urllib2 import urlopen, HTTPError, URLError
+
+import grass.script as grass
+
+class WMSBase:
+    def __init__(self, options, flags):
+        # these variables are information for destructor
+        self.temp_files_to_cleanup = []
+        self.cleanup_mask = False
+        self.cleanup_layers = False
+        
+        self.bbox     = None
+        self.temp_map = None
+        
+        if flags['c']:
+            self._getCapabilities(options)
+        else:
+            self._getMap(options, flags)  
+
+    def _getMap(self, options, flags):
+        """!Download data from WMS server and import data
+        (using GDAL library) into GRASS as a raster map."""
+
+        self._initializeParameters(options, flags)  
+
+        self.bbox     = self._computeBbox()
+        
+        self.temp_map = self._download()  
+        
+        self._createOutputMap() 
+    
+    def __del__(self):
+        # removes temporary mask, used for import transparent or warped temp_map
+        if self.cleanup_mask:
+            # clear temporary mask, which was set by module      
+            if grass.run_command('r.mask',
+                                 quiet = True,
+                                 flags = 'r') != 0:  
+                grass.fatal(_('%s failed') % 'r.mask')
+            
+            # restore original mask, if exists 
+            if grass.find_file(self.o_output + self.original_mask_suffix, element = 'cell', mapset = '.' )['name']:
+                if grass.run_command('g.copy',
+                                     quiet = True,
+                                     rast =  self.o_output + self.original_mask_suffix + ',MASK') != 0:
+                    grass.fatal(_('%s failed') % 'g.copy')
+        
+        # tries to remove temporary files, all files should be
+        # removoved before, implemented just in case of unexpected
+        # stop of module
+        for temp_file in self.temp_files_to_cleanup:
+            grass.try_remove(temp_file)
+        
+        # remove temporary created rasters
+        if self.cleanup_layers: 
+            maps = []
+            for suffix in ('.red', '.green', '.blue', '.alpha', self.original_mask_suffix):
+                rast = self.o_output + suffix
+                if grass.find_file(rast, element = 'cell', mapset = '.')['file']:
+                    maps.append(rast)
+            
+            if maps:
+                grass.run_command('g.remove',
+                                  quiet = True,
+                                  flags = 'f',
+                                  rast  = ','.join(maps))
+        
+        # deletes enviromental variable which overrides region 
+        if 'GRASS_REGION' in os.environ.keys():
+            os.environ.pop('GRASS_REGION')
+        
+    def _debug(self, fn, msg):
+        grass.debug("%s.%s: %s" %
+                    (self.__class__.__name__, fn, msg))
+        
+    def _initializeParameters(self, options, flags):
+        self._debug("_initialize_parameters", "started")
+        
+        # inicialization of module parameters (options, flags)
+        self.flags = flags 
+        if self.flags['t']:
+            self.transparent = 'TRUE'
+        else:
+            self.transparent = 'FALSE'   
+        
+        self.o_mapserver_url = options['mapserver'].strip() + "?" 
+        self.o_layers = options['layers'].strip()
+        self.o_styles = options['styles'].strip()
+        self.o_output = options['output']
+        self.o_method = options['method']
+        
+        self.o_bgcolor = options['bgcolor'].strip()
+        if self.o_bgcolor != "" and not flags["d"]:
+            grass.warning(_("Parameter bgcolor ignored, use -d flag"))
+        
+        self.o_urlparams = options['urlparams'].strip()
+        if self.o_urlparams != "" and not flags["d"]:
+            grass.warning(_("Parameter urlparams ignored, use -d flag"))
+        
+        self.o_wms_version = options['wms_version']        
+        if self.o_wms_version == "1.3.0":
+            self.projection_name = "CRS"
+        else:
+            self.projection_name = "SRS" 
+        
+        self.o_format = options['format']
+        if self.o_format == "geotiff":
+            self.mime_format = "image/geotiff"
+        elif self.o_format == "tiff":
+            self.mime_format = "image/tiff"
+        elif self.o_format == "png":
+            self.mime_format = "image/png"
+        elif self.o_format == "jpeg":
+            self.mime_format = "image/jpeg"
+            if flags['t']:
+                grass.warning(_("JPEG format does not support transparency"))
+        elif self.o_format == "gif":
+            self.mime_format = "image/gif"
+        else:
+            grass.fatal(_("Unsupported image format %s") % self.o_format)
+        
+        self.o_srs = int(options['srs'])
+        if self.o_srs <= 0:
+            grass.fatal(_("Invalid EPSG code %d") % self.o_srs)
+        
+        # read projection info
+        self.proj_location = grass.read_command('g.proj', 
+                                                flags ='jf').rstrip('\n')
+        
+        self.proj_srs = grass.read_command('g.proj', 
+                                           flags = 'jf', 
+                                           epsg = str(self.o_srs) ).rstrip('\n')
+        
+        if not self.proj_srs or not self.proj_location:
+            grass.fatal(_("Unable to get projection info"))
+        
+        # set region 
+        self.o_region = options['region']
+	if self.o_region:                 
+            if not grass.find_file(name = self.o_region, element = 'windows', mapset = '.' )['name']:
+                grass.fatal(_("Region <%s> not found") % self.o_region)
+        
+        if self.o_region:
+            s = grass.read_command('g.region',
+                                   quiet = True,
+                                   flags = 'ug',
+                                   region = self.o_region)
+            self.region = grass.parse_key_val(s, val_type = float)
+        else:
+            self.region = grass.region()
+        
+        min_tile_size = 100
+        self.o_maxcols = int(options['maxcols'])
+        if self.o_maxcols <= min_tile_size:
+            grass.fatal(_("Maxcols must be greater than 100"))
+        
+        self.o_maxrows = int(options['maxrows'])
+        if self.o_maxrows <= min_tile_size:
+            grass.fatal(_("Maxrows must be greater than 100"))
+        
+        # setting optimal tile size according to maxcols and maxrows constraint and region cols and rows      
+        self.tile_cols = int(self.region['cols'] / ceil(self.region['cols'] / float(self.o_maxcols)))
+        self.tile_rows = int(self.region['rows'] / ceil(self.region['rows'] / float(self.o_maxrows)))
+        
+        # suffix for existing mask (during overriding will be saved
+        # into raster named:self.o_output + this suffix)
+        self.original_mask_suffix = "_temp_MASK"
+        
+        # check names of temporary rasters, which module may create 
+        maps = []
+        for suffix in ('.red', '.green', '.blue', '.alpha', self.original_mask_suffix ):
+            rast = self.o_output + suffix
+            if grass.find_file(rast, element = 'cell', mapset = '.')['file']:
+                maps.append(rast)
+        
+        if len(maps) != 0:
+            grass.fatal(_("Please change output name, or change names of these rasters: %s, "
+                          "module needs to create this temporary maps during runing") % ",".join(maps))
+        
+        # default format for GDAL library
+        self.gdal_drv_format = "GTiff"
+        
+        self._debug("_initialize_parameters", "finished")
+        
+    def _getCapabilities(self, options): 
+        """!Get capabilities from WMS server
+        """
+        # download capabilities file
+        cap_url = options['mapserver'] + "service=WMS&request=GetCapabilities&version=" + options['wms_version'] 
+        try:
+            cap = urlopen(cap_url)
+        except IOError:
+            grass.fatal(_("Unable to get capabilities from '%s'") % options['mapserver'])
+        
+        cap_lines = cap.readlines()
+        for line in cap_lines: 
+            grass.message(line) 
+        
+    def _computeBbox(self):
+        """!Get region extent for WMS query (bbox)
+        """
+        self._debug("_computeBbox", "started")
+        
+        bbox = {'n' : None, 's' : None, 'e' : None, 'w' : None}
+        
+        if self.proj_srs == self.proj_location: # TODO: do it better
+            for param in bbox:
+                bbox[param] = self.region[param]
+        
+        # if location projection and wms query projection are
+        # different, corner points of region are transformed into wms
+        # projection and then bbox is created from extreme coordinates
+        # of the transformed points
+        else:
+            temp_region = self._tempfile()
+            
+            try:
+                temp_region_opened = open(temp_region, 'w')
+                temp_region_opened.write("%f %f\n%f %f\n%f %f\n%f %f\n"  %\
+                                       (self.region['e'], self.region['n'],\
+                                        self.region['w'], self.region['n'],\
+                                        self.region['w'], self.region['s'],\
+                                        self.region['e'], self.region['s'] ))
+            except IOError:
+                 grass.fatal(_("Unable to write data into tempfile"))
+            finally:           
+                temp_region_opened.close()            
+            
+            points = grass.read_command('m.proj', flags = 'd',
+                                        proj_output = self.proj_srs,
+                                        proj_input = self.proj_location,
+                                        input = temp_region) # TODO: stdin
+            grass.try_remove(temp_region)
+            if not points:
+                grass.fatal(_("Unable to determine region, %s failed") % 'm.proj')
+            
+            points = points.splitlines()
+            if len(points) != 4:
+                grass.fatal(_("Region defintion: 4 points required"))
+            
+            for point in points:
+                point = map(float, point.split("|"))
+                if not bbox['n']:
+                    bbox['n'] = point[1]
+                    bbox['s'] = point[1]
+                    bbox['e'] = point[0]
+                    bbox['w'] = point[0]
+                    continue
+                
+                if   bbox['n'] < point[1]:
+                    bbox['n'] = point[1]
+                elif bbox['s'] > point[1]:
+                    bbox['s'] = point[1]
+                
+                if   bbox['e'] < point[0]:
+                    bbox['e'] = point[0]
+                elif bbox['w'] > point[0]:
+                    bbox['w'] = point[0]  
+        
+        self._debug("_computeBbox", "finished -> %s" % bbox)
+        
+        return bbox
+
+    def _createOutputMap(self): 
+        """!Import downloaded data into GRASS, reproject data if needed
+        using gdalwarp
+        """
+        # reprojection of raster
+        if self.proj_srs != self.proj_location: # TODO: do it better
+            grass.message(_("Reprojecting raster..."))
+            temp_warpmap = self._tempfile()
+            
+            if int(os.getenv('GRASS_VERBOSE', '2')) <= 2:
+                nuldev = file(os.devnull, 'w+')
+            else:
+                nuldev = None
+            
+            # RGB rasters - alpha layer is added for cropping edges of projected raster
+            if self.temp_map_bands_num == 3:
+                ps = grass.Popen(['gdalwarp',
+                                  '-s_srs', '%s' % self.proj_srs,
+                                  '-t_srs', '%s' % self.proj_location,
+                                  '-r', self.o_method, '-dstalpha',
+                                  self.temp_map, temp_warpmap], stdout = nuldev)
+            # RGBA rasters
+            else:
+                ps = grass.Popen(['gdalwarp',
+                                  '-s_srs', '%s' % self.proj_srs,
+                                  '-t_srs', '%s' % self.proj_location,
+                                  '-r', self.o_method,
+                                  self.temp_map, temp_warpmap], stdout = nuldev)
+            ps.wait()
+            
+            if nuldev:
+                nuldev.close()
+            
+            if ps.returncode != 0:
+                grass.fatal(_('%s failed') % 'gdalwarp')
+        # raster projection is same as projection of location
+        else:
+            temp_warpmap = self.temp_map
+        
+        grass.message(_("Importing raster map into GRASS..."))
+        # importing temp_map into GRASS
+        if grass.run_command('r.in.gdal',
+                             quiet = True,
+                             input = temp_warpmap,
+                             output = self.o_output) != 0:
+            grass.fatal(_('%s failed') % 'r.in.gdal')
+        
+        # information for destructor to cleanup temp_layers, created
+        # with r.in.gdal
+        self.cleanup_layers = True
+        
+        # setting region for full extend of imported raster
+        os.environ['GRASS_REGION'] = grass.region_env(rast = self.o_output + '.red')
+        
+        # mask created from alpha layer, which describes real extend
+        # of warped layer (may not be a rectangle), also mask contains
+        # transparent parts of raster
+        if grass.find_file( self.o_output + '.alpha', element = 'cell', mapset = '.' )['name']:
+            # saving current mask (if exists) into temp raster
+            if grass.find_file('MASK', element = 'cell', mapset = '.' )['name']:
+                if grass.run_command('g.copy',
+                                     quiet = True,
+                                     rast = 'MASK,' + self.o_output + self.original_mask_suffix) != 0:    
+                    grass.fatal(_('%s failed') % 'g.copy')
+            
+            # info for destructor
+            self.cleanup_mask = True
+            if grass.run_command('r.mask',
+                                 quiet = True,
+                                 overwrite = True,
+                                 maskcats = "0",
+                                 flags = 'i',
+                                 input = self.o_output + '.alpha') != 0: 
+                grass.fatal(_('%s failed') % 'r.mask')
+        
+        if grass.run_command('r.composite',
+                             quiet = True,
+                             red = self.o_output + '.red',
+                             green = self.o_output +  '.green',
+                             blue = self.o_output + '.blue',
+                             output = self.o_output ) != 0:
+                grass.fatal(_('%s failed') % 'r.composite')
+        
+        grass.try_remove(temp_warpmap)
+        grass.try_remove(self.temp_map) 
+
+    def _tempfile(self):
+        """!Create temp_file and append list self.temp_files_to_cleanup 
+            with path of file 
+     
+        @return string path to temp_file
+        """
+        temp_file = grass.tempfile()
+        if temp_file is None:
+            grass.fatal(_("Unable to create temporary files"))
+        
+        # list of created tempfiles for destructor
+        self.temp_files_to_cleanup.append(temp_file)
+        
+        return temp_file


Property changes on: grass-addons/grass7/raster/r.in.wms2/wms_base.py
___________________________________________________________________
Added: svn:mime-type
   + text/x-python
Added: svn:eol-style
   + native

Added: grass-addons/grass7/raster/r.in.wms2/wms_drv.py
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/wms_drv.py	                        (rev 0)
+++ grass-addons/grass7/raster/r.in.wms2/wms_drv.py	2012-04-10 12:21:59 UTC (rev 51352)
@@ -0,0 +1,213 @@
+try:
+    from osgeo import gdal
+    from osgeo import gdalconst 
+except:
+    grass.fatal(_("Unable to load GDAL python bindings"))
+
+from urllib2 import urlopen
+
+import numpy as Numeric
+Numeric.arrayrange = Numeric.arange
+
+import grass.script as grass
+
+from wms_base import WMSBase
+
+class WMSDrv(WMSBase):
+    def _download(self):
+        """!Downloads data from WMS server using own driver
+        
+        @return temp_map with stored downloaded data
+        """ 
+        grass.message(_("Downloading data from WMS server..."))
+        
+        proj = self.projection_name + "=EPSG:"+ str(self.o_srs)
+        url = self.o_mapserver_url + "REQUEST=GetMap&VERSION=%s&LAYERS=%s&WIDTH=%s&HEIGHT=%s&STYLES=%s&BGCOLOR=%s&TRANSPARENT=%s" %\
+                  (self.o_wms_version, self.o_layers, self.tile_cols, self.tile_rows, self.o_styles, self.o_bgcolor, self.transparent)
+        url += "&" +proj+ "&" + "FORMAT=" + self.mime_format
+        
+        if self.o_urlparams != "":
+            url +="&" + self.o_urlparams
+        
+        cols = int(self.region['cols'])
+        rows = int(self.region['rows'])
+        
+        # computes parameters of tiles 
+        num_tiles_x = cols / self.tile_cols 
+        last_tile_x_size = cols % self.tile_cols
+        tile_x_length =  float(self.tile_cols) / float(cols ) * (self.bbox['e'] - self.bbox['w']) 
+        
+        last_tile_x = False
+        if last_tile_x_size != 0:
+            last_tile_x = True
+            num_tiles_x = num_tiles_x + 1
+        
+        num_tiles_y = rows / self.tile_rows 
+        last_tile_y_size = rows % self.tile_rows
+        tile_y_length =  float(self.tile_rows) / float(rows) * (self.bbox['n'] - self.bbox['s']) 
+        
+        last_tile_y = False
+        if last_tile_y_size != 0:
+            last_tile_y = True
+            num_tiles_y = num_tiles_y + 1
+        
+        # each tile is downloaded and written into temp_map 
+        tile_bbox = dict(self.bbox)
+        tile_bbox['e'] = self.bbox['w']  + tile_x_length
+        
+        tile_to_temp_map_size_x = self.tile_cols
+        for i_x in range(num_tiles_x):
+            # set bbox for tile i_x,i_y (E, W)
+            if i_x != 0:
+                tile_bbox['e'] += tile_x_length 
+                tile_bbox['w'] += tile_x_length            
+            
+            if i_x == num_tiles_x - 1 and last_tile_x:
+                tile_to_temp_map_size_x = last_tile_x_size 
+            
+            tile_bbox['n'] = self.bbox['n']                    
+            tile_bbox['s'] = self.bbox['n'] - tile_y_length 
+            tile_to_temp_map_size_y = self.tile_rows       
+            
+            for i_y in range(num_tiles_y):
+                # set bbox for tile i_x,i_y (N, S)
+                if i_y != 0:
+                    tile_bbox['s'] -= tile_y_length 
+                    tile_bbox['n'] -= tile_y_length
+                
+                if i_y == num_tiles_y - 1 and last_tile_y:
+                    tile_to_temp_map_size_y = last_tile_y_size 
+                
+                # bbox for tile defined
+                query_url = url + "&" + "BBOX=%s,%s,%s,%s" % (tile_bbox['w'], tile_bbox['s'], tile_bbox['e'], tile_bbox['n'])
+                
+                try: 
+                    wms_data = urlopen(query_url)
+                except IOError:
+                    grass.fatal(_("Unable to fetch data from mapserver"))
+                
+                temp_tile = self._tempfile()
+                
+                # download data into temporary file
+                try:
+                    temp_tile_opened = open(temp_tile, 'w')
+                    temp_tile_opened.write(wms_data.read())
+                except IOError:
+                    grass.fatal(_("Unable to write data into tempfile"))
+                finally:
+                    temp_tile_opened.close()
+                
+                tile_dataset_info = gdal.Open(temp_tile, gdal.GA_ReadOnly) 
+                if tile_dataset_info is None:
+                    # print error xml returned from server
+                    try:
+                        error_xml_opened = open(temp_tile, 'r')
+                        err_str = error_xml_opened.read()     
+                    except IOError:
+                        grass.fatal(_("Unable to read data from tempfile"))
+                    finally:
+                        error_xml_opened.close()
+
+                    if  err_str is not None:
+                        grass.fatal(_("WMS server error: %s") %  err_str)
+                    else:
+                        grass.fatal(_("WMS server unknown error") )
+                
+                band = tile_dataset_info.GetRasterBand(1) 
+                cell_type_func = band.__swig_getmethods__["DataType"]#??
+                bands_number_func = tile_dataset_info.__swig_getmethods__["RasterCount"]
+                
+                ##### see original r.in.wms - file gdalwarp.py line 117 ####
+                temp_tile_pct2rgb = None
+                if bands_number_func(tile_dataset_info) == 1 and band.GetRasterColorTable() is not None:
+                    # expansion of color table into bands 
+                    temp_tile_pct2rgb = self._tempfile()
+                    tile_dataset = self._pct2rgb(temp_tile, temp_tile_pct2rgb)
+                else: 
+                    tile_dataset = tile_dataset_info
+                
+                # initialization of temp_map_dataset, where all tiles are merged
+                if i_x == 0 and i_y == 0:
+                    temp_map = self._tempfile()
+                    
+                    driver = gdal.GetDriverByName(self.gdal_drv_format)
+                    metadata = driver.GetMetadata()
+                    if not metadata.has_key(gdal.DCAP_CREATE) or \
+                           metadata[gdal.DCAP_CREATE] == 'NO':
+                        grass.fatal(_('Driver %s does not supports Create() method') % drv_format)  
+                    
+                    self.temp_map_bands_num = bands_number_func(tile_dataset)
+                    temp_map_dataset = driver.Create(temp_map, int(cols), int(rows),
+                                                     self.temp_map_bands_num, cell_type_func(band));
+                
+                # tile written into temp_map
+                tile_to_temp_map = tile_dataset.ReadRaster(0, 0, tile_to_temp_map_size_x, tile_to_temp_map_size_y,
+                                                                 tile_to_temp_map_size_x, tile_to_temp_map_size_y)
+                
+                temp_map_dataset.WriteRaster(self.tile_cols * i_x, self.tile_rows * i_y,
+                                             tile_to_temp_map_size_x,  tile_to_temp_map_size_y, tile_to_temp_map) 
+                
+                tile_dataset = None
+                tile_dataset_info = None
+                grass.try_remove(temp_tile)
+                grass.try_remove(temp_tile_pct2rgb)    
+        
+        # georeferencing and setting projection of temp_map
+        projection = grass.read_command('g.proj', 
+                                        flags = 'wf',
+                                        epsg =self.o_srs).rstrip('\n')
+        temp_map_dataset.SetProjection(projection)
+        
+        pixel_x_length = (self.bbox['e'] - self.bbox['w']) / int(cols)
+        pixel_y_length = (self.bbox['s'] - self.bbox['n']) / int(rows)
+        geo_transform = [ self.bbox['w'] , pixel_x_length  , 0.0 ,  self.bbox['n'] , 0.0 , pixel_y_length ] 
+        temp_map_dataset.SetGeoTransform(geo_transform )
+        temp_map_dataset = None
+        
+        return temp_map
+    
+    def _pct2rgb(self, src_filename, dst_filename):
+        """!Create new dataset with data in dst_filename with bands according to src_filename 
+        raster color table - modified code from gdal utility pct2rgb
+        
+        @return new dataset
+        """  
+        out_bands = 4
+        band_number = 1
+        
+        # open source file
+        src_ds = gdal.Open(src_filename)
+        if src_ds is None:
+            grass.fatal(_('Unable to open %s ' % src_filename))
+            
+        src_band = src_ds.GetRasterBand(band_number)
+        
+        # Build color table
+        lookup = [ Numeric.arrayrange(256), 
+                   Numeric.arrayrange(256), 
+                   Numeric.arrayrange(256), 
+                   Numeric.ones(256)*255 ]
+        
+        ct = src_band.GetRasterColorTable()	
+        if ct is not None:
+            for i in range(min(256,ct.GetCount())):
+                entry = ct.GetColorEntry(i)
+                for c in range(4):
+                    lookup[c][i] = entry[c]
+        
+        # create the working file
+        gtiff_driver = gdal.GetDriverByName(self.gdal_drv_format)
+        tif_ds = gtiff_driver.Create(dst_filename,
+                                     src_ds.RasterXSize, src_ds.RasterYSize, out_bands)
+        
+        # do the processing one scanline at a time
+        for iY in range(src_ds.RasterYSize):
+            src_data = src_band.ReadAsArray(0,iY,src_ds.RasterXSize,1)
+            
+            for iBand in range(out_bands):
+                band_lookup = lookup[iBand]
+                
+                dst_data = Numeric.take(band_lookup,src_data)
+                tif_ds.GetRasterBand(iBand+1).WriteArray(dst_data,0,iY)
+        
+        return tif_ds       


Property changes on: grass-addons/grass7/raster/r.in.wms2/wms_drv.py
___________________________________________________________________
Added: svn:mime-type
   + text/x-python
Added: svn:eol-style
   + native

Added: grass-addons/grass7/raster/r.in.wms2/wms_gdal_drv.py
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/wms_gdal_drv.py	                        (rev 0)
+++ grass-addons/grass7/raster/r.in.wms2/wms_gdal_drv.py	2012-04-10 12:21:59 UTC (rev 51352)
@@ -0,0 +1,131 @@
+import os
+
+try:
+    from osgeo import gdal
+    from osgeo import gdalconst
+except:
+    grass.fatal(_("Unable to load GDAL Python bindings"))
+
+import xml.etree.ElementTree as etree
+
+import grass.script as grass
+
+from wms_base import WMSBase
+
+class NullDevice():
+    def write(self, s):
+        pass
+
+class WMSGdalDrv(WMSBase):
+    def _createXML(self):
+        """!Create XML for GDAL WMS driver
+        
+        @return path to XML file
+        """
+        self._debug("_createXML", "started")
+        
+        gdal_wms = etree.Element("GDAL_WMS")
+        service = etree.SubElement(gdal_wms, "Service")
+        name = etree.Element("name")
+        service.set("name","WMS")
+        
+        version = etree.SubElement(service, "Version")
+        version.text =self.o_wms_version
+        
+        server_url = etree.SubElement(service, "ServerUrl")
+        server_url.text =self.o_mapserver_url
+        
+        srs = etree.SubElement(service, self.projection_name)   
+        srs.text = 'EPSG:' + str(self.o_srs)
+        
+        image_format = etree.SubElement(service, "ImageFormat")
+        image_format.text = self.mime_format
+        
+        image_format = etree.SubElement(service, "Transparent")
+        image_format.text = self.transparent
+        
+        layers = etree.SubElement(service, "Layers")
+        layers.text = self.o_layers
+        
+        styles = etree.SubElement(service, "Styles")
+        styles.text = self.o_styles
+        
+        data_window = etree.SubElement(gdal_wms, "DataWindow")
+        
+        upper_left_x = etree.SubElement(data_window, "UpperLeftX")
+        upper_left_x.text = str(self.bbox['w']) 
+        
+        upper_left_y = etree.SubElement(data_window, "UpperLeftY")
+        upper_left_y.text = str(self.bbox['n']) 
+        
+        lower_right_x = etree.SubElement(data_window, "LowerRightX")
+        lower_right_x.text = str(self.bbox['e']) 
+        
+        lower_right_y = etree.SubElement(data_window, "LowerRightY")
+        lower_right_y.text = str(self.bbox['s'])
+        
+        size_x = etree.SubElement(data_window, "SizeX")
+        size_x.text = str(self.region['cols']) 
+        
+        size_y = etree.SubElement(data_window, "SizeY")
+        size_y.text = str(self.region['rows']) 
+        
+        # RGB + alpha
+        self.temp_map_bands_num = 4
+        block_size_x = etree.SubElement(gdal_wms, "BandsCount")
+        block_size_x.text = str(self.temp_map_bands_num)
+        
+        block_size_x = etree.SubElement(gdal_wms, "BlockSizeX")
+        block_size_x.text = str(self.tile_cols) 
+        
+        block_size_y = etree.SubElement(gdal_wms, "BlockSizeY")
+        block_size_y.text = str(self.tile_rows)
+        
+        xml_file = self._tempfile()
+        
+        etree.ElementTree(gdal_wms).write(xml_file)
+        
+        self._debug("_createXML", "finished -> %s" % xml_file)
+        
+        return xml_file
+    
+    def _download(self):
+        """!Downloads data from WMS server using GDAL WMS driver
+        
+        @return temp_map with stored downloaded data
+        """
+        grass.message("Downloading data from WMS server...")
+        self._debug("_download", "started")
+        
+        temp_map = self._tempfile()
+        
+        xml_file = self._createXML()
+        wms_dataset = gdal.Open(xml_file, gdal.GA_ReadOnly)
+        grass.try_remove(xml_file)
+        if wms_dataset is None:
+            grass.fatal(_("Unable to open GDAL WMS driver"))
+        
+        self._debug("_download", "GDAL dataset created")
+        
+        driver = gdal.GetDriverByName(self.gdal_drv_format)
+        if driver is None:
+            grass.fatal(_("Unable to find %s driver" % format))
+        
+        metadata = driver.GetMetadata()
+        if not metadata.has_key(gdal.DCAP_CREATECOPY) or \
+           metadata[gdal.DCAP_CREATECOPY] == 'NO':
+            grass.fatal(_('Driver %s supports CreateCopy() method.') % self.gdal_drv_name)
+        
+        self._debug("_download", "calling GDAL CreateCopy...")
+        
+        temp_map_dataset = driver.CreateCopy(temp_map, wms_dataset, 0)
+        
+        if temp_map_dataset is None:
+            grass.fatal(_("Incorrect WMS query"))
+        
+        temp_map_dataset  = None
+        wms_dataset = None
+        
+        self._debug("_download", "finished")
+        
+        return temp_map


Property changes on: grass-addons/grass7/raster/r.in.wms2/wms_gdal_drv.py
___________________________________________________________________
Added: svn:mime-type
   + text/x-python
Added: svn:eol-style
   + native



More information about the grass-commit mailing list