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

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Oct 9 15:02:55 PDT 2012


Author: turek
Date: 2012-10-09 15:02:55 -0700 (Tue, 09 Oct 2012)
New Revision: 53348

Modified:
   grass-addons/grass7/raster/r.in.wms2/r.in.wms2.html
   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:
r.in.wms2: WMTS support added

Modified: grass-addons/grass7/raster/r.in.wms2/r.in.wms2.html
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/r.in.wms2.html	2012-10-09 21:50:21 UTC (rev 53347)
+++ grass-addons/grass7/raster/r.in.wms2/r.in.wms2.html	2012-10-09 22:02:55 UTC (rev 53348)
@@ -1,8 +1,9 @@
 <h2>DESCRIPTION</h2>
 
-<em>r.in.wms</em> handles all of downloading and importing raster
+<em>r.in.wms2</em> handles all of downloading and importing raster
 data from an <a href="http://www.opengeospatial.org/standards/wms">OGC
-WMS</a> web mapping server. It need only be told the desired data to
+WMS</a> and <a href="http://www.opengeospatial.org/standards/wmts">OGC
+WMTS</a> web mapping servers. It need only be told the desired data to
 collect (bounds and resolution) via a region, the server to get the
 data from, and the layer or layers to get. It downloads the data in
 tiles, reprojects it, imports it, and patches it back together.
@@ -25,39 +26,43 @@
 <h3>General Get Capabilities Request</h3>
 
 <div class="code"><pre>
-r.in.wms -c mapserver=http://wms.cuzk.cz/wms.asp
+r.in.wms2 -c url=http://wms.cuzk.cz/wms.asp
 </pre></div>
 
-<h3>CUZK download</h3>
+<h3>>Download raster data from WMS server (GetMap request)</h3>
 
 
 World extend data:
 
 <div class="code"><pre>
-r.in.wms mapserver=http://iceds.ge.ucl.ac.uk/cgi-bin/icedswms layers=bluemarble,landsat_1_01 styles=default,default output=landsat srs=4326 format=png 
+r.in.wms2 url=http://iceds.ge.ucl.ac.uk/cgi-bin/icedswms layers=bluemarble,landsat_1_01 styles=default,default output=landsat srs=4326 format=png 
 </pre></div>
 * Server supports only WMS 1.1.1 <br>
 
 <div class="code"><pre>
-r.in.wms mapserver=http://132.156.97.59/cgi-bin/worldmin_en-ca_ows layers=GSC:WORLD_PrecambrianDomains output=pokus srs=4326 format=jpeg 
+r.in.wms2 url=http://132.156.97.59/cgi-bin/worldmin_en-ca_ows layers=GSC:WORLD_PrecambrianDomains output=pokus srs=4326 format=jpeg 
 </pre></div>
 * Server supports only WMS 1.1.1 
 <br>
 <br>
 
+<div class="code"><pre>
+r.in.wms2 url=http://gpp3-wxs.ign.fr/yourAPIkey/geoportail/wmts layers=ORTHOIMAGERY.ORTHOPHOTOS output=orthophoto srs=3857 format=jpeg driver=WMTS_GRASS style=normal password=* username=*
+</pre></div>
+* Username, password and API key can be get at <a href="http://api.ign.fr/">IGN API</a> website
+<br>
+<br>
+
 Data in extend of Czech Republic:
 
 <div class="code"><pre>
-r.in.wms output=kn mapserver=http://wms.cuzk.cz/wms.asp layers=prehledka_kraju-linie srs=4326 format=png
+r.in.wms2 output=kn url=http://wms.cuzk.cz/wms.asp layers=prehledka_kraju-linie srs=4326 format=png
 </pre></div>
 
+<div class="code"><pre>
+r.in.wms2 url=http://geoportal.cuzk.cz/WMTS_ORTOFOTO/WMTService.aspx layers=orto output=ortofoto srs=3857 format=jpeg driver=WMTS_GRASS style=default
+</pre></div>
 
-<h2>TODO</h2>
-
-<ul>
-  <li>Implement Tiled WMS</li>
-</ul>
-
 <h2>SEE ALSO</h2>
 
 <em>

Modified: grass-addons/grass7/raster/r.in.wms2/r.in.wms2.py
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/r.in.wms2.py	2012-10-09 21:50:21 UTC (rev 53347)
+++ grass-addons/grass7/raster/r.in.wms2/r.in.wms2.py	2012-10-09 22:02:55 UTC (rev 53348)
@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 """
-MODULE:    r.in.wms
+MODULE:    r.in.wms2
 
 AUTHOR(S): Stepan Turek <stepan.turek AT seznam.cz>
 
@@ -20,7 +20,7 @@
 #%end
 
 #%option
-#% key: mapserver
+#% key: url
 #% type: string
 #% description:URL of WMS server 
 #% required: yes
@@ -42,6 +42,7 @@
 #% key: srs
 #% type: integer
 #% description: EPSG number of source projection for request 
+#% answer:4326 
 #% guisection: Request properties
 #%end
 
@@ -98,11 +99,25 @@
 #%option
 #% key: urlparams
 #% type:string
-#% description: Addition query parameters for server (only with 'd' flag)
+#% description: Additional query parameters for server
 #% guisection: Request properties
 #%end
 
 #%option
+#% key: username
+#% type:string
+#% description: Username for server connection
+#% guisection: Request properties
+#%end
+
+#%option
+#% key: password
+#% type:string
+#% description: Password for server connection
+#% guisection: Request properties
+#%end
+
+#%option
 #% key: styles
 #% type: string
 #% description: Styles to request from map server
@@ -113,7 +128,7 @@
 #%option
 #% key: bgcolor
 #% type: string
-#% description: Color of map background (only with 'd' flag)
+#% description: Color of map background
 #% guisection: Map style
 #%end
 
@@ -130,23 +145,27 @@
 #% suppress_required: yes
 #%end
 
-#%flag
-#% key: d
-#% description: Do not use GDAL WMS driver
+#%option
+#% key: driver
+#% type:string
+#% description: Driver for communication with server
+#% options:WMS_GDAL, WMS_GRASS, WMTS_GRASS
+#% answer:WMS_GDAL
 #%end
 
+
 import os
 import sys
-sys.path.insert(1, os.path.join(os.path.dirname(sys.path[0]), 'etc', 'r.in.wms'))
+sys.path.insert(1, os.path.join(os.path.dirname(sys.path[0]), 'etc', 'r.in.wms2'))
 
 import grass.script as grass
 
 def main():
-    if flags['d']:
-        grass.debug("Using own driver")
+    if 'GRASS' in options['driver']:
+        grass.debug("Using GRASS driver")
         from wms_drv import WMSDrv
         wms = WMSDrv()
-    else:
+    elif 'GDAL' in options['driver']:
         grass.debug("Using GDAL WMS driver")
         from wms_gdal_drv import WMSGdalDrv
         wms = WMSGdalDrv()

Modified: grass-addons/grass7/raster/r.in.wms2/wms_base.py
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/wms_base.py	2012-10-09 21:50:21 UTC (rev 53347)
+++ grass-addons/grass7/raster/r.in.wms2/wms_base.py	2012-10-09 22:02:55 UTC (rev 53348)
@@ -1,8 +1,8 @@
 import os
 from   math import ceil
 
-import xml.etree.ElementTree as etree
 from urllib2 import urlopen, HTTPError, URLError
+from httplib import HTTPException
 
 import grass.script as grass
 
@@ -13,7 +13,9 @@
         self.cleanup_mask   = False
         self.cleanup_layers = False
         
-        self.bbox     = None
+        self.params = {}
+        self.tile_size = {'bbox' : None}
+
         self.temp_map = None
         
     def __del__(self):
@@ -26,14 +28,14 @@
                 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.find_file(self.params['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:
+                                     rast =  self.params['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
+        # removed before, implemented just in case of unexpected
         # stop of module
         for temp_file in self.temp_files_to_cleanup:
             grass.try_remove(temp_file)
@@ -42,7 +44,7 @@
         if self.cleanup_layers: 
             maps = []
             for suffix in ('.red', '.green', '.blue', '.alpha', self.original_mask_suffix):
-                rast = self.o_output + suffix
+                rast = self.params['output'] + suffix
                 if grass.find_file(rast, element = 'cell', mapset = '.')['file']:
                     maps.append(rast)
             
@@ -52,7 +54,7 @@
                                   flags = 'f',
                                   rast  = ','.join(maps))
         
-        # deletes enviromental variable which overrides region 
+        # deletes environmental variable which overrides region 
         if 'GRASS_REGION' in os.environ.keys():
             os.environ.pop('GRASS_REGION')
         
@@ -63,100 +65,109 @@
     def _initializeParameters(self, options, flags):
         self._debug("_initialize_parameters", "started")
         
-        # inicialization of module parameters (options, flags)
-        self.flags = flags 
-        if self.flags['o']:
-            self.transparent = 'FALSE'
+        # initialization of module parameters (options, flags)
+
+        self.params['driver'] = options['driver']
+
+        self.flags = flags
+
+        if self.flags['o'] and 'WMS' not in self.params['driver']:
+            grass.warning(_("Flag '%s' is relevant only for WMS.") % 'o')
+        elif self.flags['o']:
+            self.params['transparent'] = 'FALSE'
         else:
-            self.transparent = 'TRUE'   
+            self.params['transparent'] = 'TRUE'   
+
+        for key in ['url', 'layers', 'styles', 'output', 'method']:
+            self.params[key] = options[key].strip()
+
+        for key in ['password', 'username', 'urlparams']:
+            self.params[key] = options[key] 
+            if self.params[key] != "" and 'GRASS' not in self.params['driver']:
+                grass.warning(_("Parameter '%s' is relevant only for %s drivers.") % (key, '*_GRASS'))
         
-        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"
+        self.params['bgcolor'] = options['bgcolor'].strip()
+        if self.params['bgcolor'] != "" and 'WMS_GRASS' not in self.params['driver']:
+            grass.warning(_("Parameter '%s' is relevant only for %s driver.") % ('bgcolor', 'WMS_GRASS'))
+                
+        self.params['wms_version'] = options['wms_version']        
+        if self.params['wms_version'] == "1.3.0":
+            self.params['proj_name'] = "CRS"
         else:
-            self.projection_name = "SRS" 
+            self.params['proj_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  options['format'] == "geotiff":
+            self.params['format'] = "image/geotiff"
+        elif options['format'] == "tiff":
+            self.params['format'] = "image/tiff"
+        elif options['format'] == "png":
+            self.params['format'] = "image/png"
+        elif  options['format'] == "jpeg":
+            self.params['format'] = "image/jpeg"
             if flags['o']:
                 grass.warning(_("JPEG format does not support transparency"))
-        elif self.o_format == "gif":
-            self.mime_format = "image/gif"
+        elif self.params['format'] == "gif":
+            self.params['format'] = "image/gif"
         else:
-            self.mime_format = self.o_format
+            self.params['format'] = self.params['format']
         
-        self.o_srs = int(options['srs'])
-        if self.o_srs <= 0:
-            grass.fatal(_("Invalid EPSG code %d") % self.o_srs)
+        self.params['srs'] = int(options['srs'])
+        if self.params['srs'] <= 0:
+            grass.fatal(_("Invalid EPSG code %d") % self.params['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 self.params['srs'] in [3857, 900913]:
+            # HACK: epsg 3857 def: http://spatialreference.org/ref/sr-org/7483/
+            # g.proj can return: ...+a=6378137 +rf=298.257223563... (WGS84 elipsoid def instead of sphere), it can make 20km shift in Y, when raster is transformed
+            # needed to be tested on more servers
+            self.proj_srs = '+proj=merc +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +no_defs +a=6378137 +b=6378137 +nadgrids=@null +to_meter=1 +wktext'
+        else:
+            self.proj_srs = grass.read_command('g.proj', 
+                                               flags = 'jf', 
+                                               epsg = str(self.params['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)
+        self.params['region'] = options['region']
+        if self.params['region']:                 
+            if not grass.find_file(name = self.params['region'], element = 'windows', mapset = '.' )['name']:
+                grass.fatal(_("Region <%s> not found") % self.params['region'])
         
-        if self.o_region:
+        if self.params['region']:
             s = grass.read_command('g.region',
                                    quiet = True,
                                    flags = 'ug',
-                                   region = self.o_region)
+                                   region = self.params['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:
+        maxcols = int(options['maxcols'])
+        if 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:
+        maxrows = int(options['maxrows'])
+        if 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)))
+        self.tile_size['cols'] = int(self.region['cols'] / ceil(self.region['cols'] / float(maxcols)))
+        self.tile_size['rows'] = int(self.region['rows'] / ceil(self.region['rows'] / float(maxrows)))
         
         # suffix for existing mask (during overriding will be saved
-        # into raster named:self.o_output + this suffix)
+        # into raster named:self.params['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
+            rast = self.params['output'] + suffix
             if grass.find_file(rast, element = 'cell', mapset = '.')['file']:
                 maps.append(rast)
         
@@ -177,20 +188,35 @@
 
         self.bbox     = self._computeBbox()
         
-        self.temp_map = self._download()  
-        
+        self.temp_map = self._download()
+
+        if not self.temp_map:
+            return
+
         self._createOutputMap() 
     
-    def GetCapabilities(self, options): 
-        """!Get capabilities from WMS server
+    def _fetchCapabilities(self, options): 
+        """!Download capabilities from WMS server
         """
         # download capabilities file
-        cap_url = options['mapserver'] + "?service=WMS&request=GetCapabilities&version=" + options['wms_version'] 
+        cap_url = options['url']
+
+        if 'WMTS' in options['driver']:
+            cap_url += "?SERVICE=WMTS&REQUEST=GetCapabilities&VERSION=1.0.0"
+        else:
+            cap_url += "?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'])
+        except (IOError, HTTPError):
+            grass.fatal(_("Unable to get capabilities from '%s'") % options['url'])
         
+        return cap
+
+    def GetCapabilities(self, options): 
+        """!Get capabilities from WMS server
+        """
+        cap  = self._fetchCapabilities(options)
         cap_lines = cap.readlines()
         for line in cap_lines: 
             print line 
@@ -228,7 +254,7 @@
                  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,
@@ -239,10 +265,13 @@
             
             points = points.splitlines()
             if len(points) != 4:
-                grass.fatal(_("Region defintion: 4 points required"))
+                grass.fatal(_("Region definition: 4 points required"))
             
             for point in points:
-                point = map(float, point.split("|"))
+                try:
+                    point = map(float, point.split("|"))
+                except ValueError:
+                    grass.fatal(_('Reprojection of region using m.proj failed.'))
                 if not bbox['maxy']:
                     bbox['maxy'] = point[1]
                     bbox['miny'] = point[1]
@@ -263,16 +292,9 @@
         self._debug("_computeBbox", "finished -> %s" % bbox)
 
         # Ordering of coordinates axis of geographic coordinate
-        # systems in WMS 1.3.0 is fliped. If self.flip_coords is 
+        # systems in WMS 1.3.0 is flipped. If  self.tile_size['flip_coords'] is 
         # True, coords in bbox need to be flipped in WMS query.
 
-        self.flip_coords = False  
-        hasLongLat = self.proj_srs.find("+proj=longlat")   
-        hasLatLong = self.proj_srs.find("+proj=latlong")   
-
-        if (hasLongLat != -1 or hasLatLong != -1) and self.o_wms_version == "1.3.0":
-            self.flip_coords = True
-
         return bbox
 
     def _createOutputMap(self): 
@@ -289,19 +311,20 @@
             else:
                 nuldev = None
             
+            #"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
             # 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',
+                                  '-r', self.params['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,
+                                  '-r', self.params['method'],
                                   self.temp_map, temp_warpmap], stdout = nuldev)
             ps.wait()
             
@@ -319,7 +342,7 @@
         if grass.run_command('r.in.gdal',
                              quiet = True,
                              input = temp_warpmap,
-                             output = self.o_output) != 0:
+                             output = self.params['output']) != 0:
             grass.fatal(_('%s failed') % 'r.in.gdal')
         
         # information for destructor to cleanup temp_layers, created
@@ -327,17 +350,17 @@
         self.cleanup_layers = True
         
         # setting region for full extend of imported raster
-        os.environ['GRASS_REGION'] = grass.region_env(rast = self.o_output + '.red')
+        os.environ['GRASS_REGION'] = grass.region_env(rast = self.params['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']:
+        if grass.find_file( self.params['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:    
+                                     rast = 'MASK,' + self.params['output'] + self.original_mask_suffix) != 0:    
                     grass.fatal(_('%s failed') % 'g.copy')
             
             # info for destructor
@@ -347,39 +370,20 @@
                                  overwrite = True,
                                  maskcats = "0",
                                  flags = 'i',
-                                 input = self.o_output + '.alpha') != 0: 
+                                 input = self.params['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:
+                             red = self.params['output'] + '.red',
+                             green = self.params['output'] +  '.green',
+                             blue = self.params['output'] + '.blue',
+                             output = self.params['output'] ) != 0:
                 grass.fatal(_('%s failed') % 'r.composite')
         
         grass.try_remove(temp_warpmap)
         grass.try_remove(self.temp_map) 
 
-    def _flipBbox(self, bbox):
-        """ 
-        flips items in dictionary 
-        value flips between this keys:
-        maxy -> maxx
-        maxx -> maxy
-        miny -> minx
-        minx -> miny
-        @return copy of bbox with fliped cordinates
-        """  
-        temp_bbox = dict(bbox)
-        new_bbox = {}
-        new_bbox['maxy'] = temp_bbox['maxx']
-        new_bbox['miny'] = temp_bbox['minx']
-        new_bbox['maxx'] = temp_bbox['maxy']
-        new_bbox['minx'] = temp_bbox['miny']
-
-        return new_bbox
-
     def _tempfile(self):
         """!Create temp_file and append list self.temp_files_to_cleanup 
             with path of file 

Modified: grass-addons/grass7/raster/r.in.wms2/wms_drv.py
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/wms_drv.py	2012-10-09 21:50:21 UTC (rev 53347)
+++ grass-addons/grass7/raster/r.in.wms2/wms_drv.py	2012-10-09 22:02:55 UTC (rev 53348)
@@ -1,3 +1,8 @@
+
+import grass.script as grass 
+from urllib2 import urlopen, HTTPError
+from httplib import HTTPException
+
 try:
     from osgeo import gdal
     from osgeo import gdalconst 
@@ -4,172 +9,146 @@
 except:
     grass.fatal(_("Unable to load GDAL python bindings"))
 
-from urllib2 import urlopen
+import urllib2
+import xml.etree.ElementTree as etree
 
 import numpy as Numeric
 Numeric.arrayrange = Numeric.arange
 
-import grass.script as grass
-
+from math import pi, floor
 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
+        @return temp_map with 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['maxx'] - self.bbox['minx']) 
-        
-        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['maxy'] - self.bbox['miny']) 
-        
-        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['maxx'] = self.bbox['minx']  + 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['maxx'] += tile_x_length 
-                tile_bbox['minx'] += 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['maxy'] = self.bbox['maxy']                    
-            tile_bbox['miny'] = self.bbox['maxy'] - 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['miny'] -= tile_y_length 
-                    tile_bbox['maxy'] -= tile_y_length
-                
-                if i_y == num_tiles_y - 1 and last_tile_y:
-                    tile_to_temp_map_size_y = last_tile_y_size 
-                
-                if self.flip_coords:
-                    # flips coordinates if WMS strandard is 1.3.0 and 
-                    # projection is geographic (see:wms_base.py _computeBbox)
-                    query_bbox = dict(self._flipBbox(tile_bbox))
+        # initialize correct manager according to chosen OGC service
+        if self.params['driver'] == 'WMTS_GRASS':
+            req_mgr = WMTSRequestMgr(self.params, self.bbox, self.region, self.proj_srs, self._fetchCapabilities(self.params))
+        elif self.params['driver'] == 'WMS_GRASS':
+            req_mgr = WMSRequestMgr(self.params, self.bbox, self.region, self.tile_size, self.proj_srs)
+
+        # set password and username to urlib2
+        if self.params['username']:
+            passman = urllib2.HTTPPasswordMgrWithDefaultRealm()
+            passman.add_password(None, self.params['url'], self.params['username'], self.params['password'])
+            auth_handler = urllib2.HTTPBasicAuthHandler(passman)
+            opener = urllib2.build_opener(auth_handler)
+            urllib2.install_opener(opener)
+
+        # get information about size in pixels and bounding box of raster, where all tiles will be joined
+        map_region = req_mgr.GetMapRegion()
+
+        init = True
+        temp_map = None
+
+        # iterate through all tiles and download them
+        while True:
+
+            # get url for request the tile and information for placing the tile into raster with other tiles
+            tile = req_mgr.GetNextTile()
+
+            # if last tile has been already downloaded 
+            if not tile:
+                break
+
+            # url for request the tile
+            query_url = tile[0]
+
+            # the tile size and offset in pixels for placing it into raster where tiles are joined
+            tileRef = tile[1]
+
+            grass.debug(query_url)
+            try: 
+                wms_data = urllib2.urlopen(query_url)
+            except (IOError, HTTPException), e:
+                if HTTPError == type(e) and e.code == 401:
+                    grass.fatal(_("Authorization failed"))           
                 else:
-                    query_bbox = tile_bbox
+                    grass.fatal(_("Unable to fetch data from server")) 
 
-                query_url = url + "&" + "BBOX=%s,%s,%s,%s" % ( query_bbox['minx'],  query_bbox['miny'],  query_bbox['maxx'],  query_bbox['maxy'])
-                grass.debug(query_url)
-                try: 
-                    wms_data = urlopen(query_url)
-                except IOError:
-                    grass.fatal(_("Unable to fetch data from mapserver"))
+            temp_tile = self._tempfile()
                 
-                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()
                 
-                # download data into temporary file
+            tile_dataset_info = gdal.Open(temp_tile, gdal.GA_ReadOnly) 
+            if tile_dataset_info is None:
+                # print error xml returned from server
                 try:
-                    temp_tile_opened = open(temp_tile, 'w')
-                    temp_tile_opened.write(wms_data.read())
+                    error_xml_opened = open(temp_tile, 'r')
+                    err_str = error_xml_opened.read()     
                 except IOError:
-                    grass.fatal(_("Unable to write data into tempfile"))
+                    grass.fatal(_("Unable to read data from 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()
+                    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") )
+                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"]
+            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
+            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()
+            # initialization of temp_map_dataset, where all tiles are merged
+            if init:
+                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)  
+                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));
+                self.temp_map_bands_num = bands_number_func(tile_dataset)
+                temp_map_dataset = driver.Create(temp_map, map_region['cols'], map_region['rows'],
+                                                 self.temp_map_bands_num, cell_type_func(band))
+                init = False
                 
-                # 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)
+            # tile is written into temp_map
+            tile_to_temp_map = tile_dataset.ReadRaster(0, 0, tileRef['sizeX'], tileRef['sizeY'],
+                                                             tileRef['sizeX'], tileRef['sizeY'])
                 
-                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) 
+            temp_map_dataset.WriteRaster(tileRef['t_cols_offset'], tileRef['t_rows_offset'],
+                                         tileRef['sizeX'],  tileRef['sizeY'], tile_to_temp_map) 
                 
-                tile_dataset = None
-                tile_dataset_info = None
-                grass.try_remove(temp_tile)
-                grass.try_remove(temp_tile_pct2rgb)    
+            tile_dataset = None
+            tile_dataset_info = None
+            grass.try_remove(temp_tile)
+            grass.try_remove(temp_tile_pct2rgb)    
         
+        if not temp_map:
+            return temp_map
         # georeferencing and setting projection of temp_map
         projection = grass.read_command('g.proj', 
                                         flags = 'wf',
-                                        epsg =self.o_srs).rstrip('\n')
+                                        epsg =self.params['srs']).rstrip('\n')
         temp_map_dataset.SetProjection(projection)
-        
 
+        pixel_x_length = (map_region['maxx'] - map_region['minx']) / int(map_region['cols'])
+        pixel_y_length = (map_region['miny'] - map_region['maxy']) / int(map_region['rows'])
 
-        pixel_x_length = (self.bbox['maxx'] - self.bbox['minx']) / int(cols)
-        pixel_y_length = (self.bbox['miny'] - self.bbox['maxy']) / int(rows)
-        geo_transform = [ self.bbox['minx'] , pixel_x_length  , 0.0 ,  self.bbox['maxy'] , 0.0 , pixel_y_length ] 
+        geo_transform = [map_region['minx'] , pixel_x_length  , 0.0 , map_region['maxy'] , 0.0 , pixel_y_length] 
         temp_map_dataset.SetGeoTransform(geo_transform )
         temp_map_dataset = None
         
@@ -220,3 +199,405 @@
                 tif_ds.GetRasterBand(iBand+1).WriteArray(dst_data,0,iY)
         
         return tif_ds       
+
+class BaseRequestMgr:
+    """!Base class for request managers. 
+    """
+    def _isGeoProj(self, proj):
+        """!Is it geographic projection? 
+        """       
+        if (proj.find("+proj=latlong")  != -1 or \
+            proj.find("+proj=longlat")  != -1):
+
+            return True
+        return False
+
+    def find(self, etreeElement, tag, ns = None):
+        """!Wraper for etree element find method. 
+        """
+        if not ns:
+            res = etreeElement.find(tag)
+        else:
+            res = etreeElement.find(ns(tag))
+
+        if res is None:
+            grass.fatal(_("Unable to parse capabilities file. \n Tag '%s' was not found.") % tag)
+
+        return res
+
+    def findall(self, etreeElement, tag, ns = None):
+        """!Wraper for etree element findall method. 
+        """
+        if not ns:
+            res = etreeElement.findall(tag)
+        else:
+            res = etreeElement.findall(ns(tag))
+
+        if not res:
+            grass.fatal(_("Unable to parse capabilities file. \n Tag '%s' was not found.") % tag)
+
+        return res
+
+class WMSRequestMgr(BaseRequestMgr):
+    def __init__(self, params, bbox, region, tile_size, proj_srs, cap_file = None):
+        """!Initialize data needed for iteration through tiles.
+        """
+        self.version = params['wms_version']
+        proj = params['proj_name'] + "=EPSG:"+ str(params['srs'])
+        self.url = params['url'] + ("?SERVICE=WMS&REQUEST=GetMap&VERSION=%s&LAYERS=%s&WIDTH=%s&HEIGHT=%s&STYLES=%s&BGCOLOR=%s&TRANSPARENT=%s" % \
+                  (params['wms_version'], params['layers'], tile_size['cols'], tile_size['rows'], params['styles'], \
+                   params['bgcolor'], params['transparent']))
+        self.url += "&" +proj+ "&" + "FORMAT=" + params['format']
+        
+        self.bbox = bbox
+        self.proj_srs = proj_srs
+        self.tile_rows = tile_size['rows']
+        self.tile_cols = tile_size['cols']
+
+        if params['urlparams'] != "":
+            self.url += "&" + params['urlparams']
+        
+        cols = int(region['cols'])
+        rows = int(region['rows'])
+        
+        # computes parameters of tiles 
+        self.num_tiles_x = cols / self.tile_cols 
+        self.last_tile_x_size = cols % self.tile_cols
+        self.tile_length_x =  float(self.tile_cols) / float(cols) * (self.bbox['maxx'] - self.bbox['minx']) 
+        
+        self.last_tile_x = False
+        if self.last_tile_x_size != 0:
+            self.last_tile_x = True
+            self.num_tiles_x = self.num_tiles_x + 1
+        
+        self.num_tiles_y = rows / self.tile_rows 
+        self.last_tile_y_size = rows % self.tile_rows
+        self.tile_length_y =  float(self.tile_rows) / float(rows) * (self.bbox['maxy'] - self.bbox['miny']) 
+        
+        self.last_tile_y = False
+        if self.last_tile_y_size != 0:
+            self.last_tile_y = True
+            self.num_tiles_y = self.num_tiles_y + 1
+        
+        self.tile_bbox = dict(self.bbox)
+        self.tile_bbox['maxx'] = self.bbox['minx']  + self.tile_length_x
+
+        self.i_x = 0 
+        self.i_y = 0
+
+        self.map_region = self.bbox
+        self.map_region['cols'] = cols
+        self.map_region['rows'] = rows
+
+    def GetMapRegion(self):
+        """!Get size in pixels and bounding box of raster where all tiles will be merged.
+        """
+        return self.map_region
+
+    def GetNextTile(self):
+        """!Get url for request the tile from server and information for merging the tile with other tiles
+        """
+
+        tileRef = {}
+
+        if self.i_x >= self.num_tiles_x:
+            return None
+            
+        tileRef['sizeX'] = self.tile_cols
+        if self.i_x == self.num_tiles_x - 1 and self.last_tile_x:
+            tileRef['sizeX'] = self.last_tile_x_size
+         
+        # set bbox for tile (N, S)
+        if self.i_y != 0:
+            self.tile_bbox['miny'] -= self.tile_length_y 
+            self.tile_bbox['maxy'] -= self.tile_length_y
+        else:
+            self.tile_bbox['maxy'] = self.bbox['maxy']
+            self.tile_bbox['miny'] = self.bbox['maxy'] - self.tile_length_y
+
+        tileRef['sizeY'] = self.tile_rows
+        if self.i_y == self.num_tiles_y - 1 and self.last_tile_y:
+            tileRef['sizeY'] = self.last_tile_y_size 
+
+        if self._isGeoProj(self.proj_srs) and self.version == "1.3.0":                
+            query_bbox = self._flipBbox(self.tile_bbox, self.proj_srs)
+        else:
+            query_bbox = self.tile_bbox
+        query_url = self.url + "&" + "BBOX=%s,%s,%s,%s" % ( query_bbox['minx'],  query_bbox['miny'],  query_bbox['maxx'],  query_bbox['maxy'])
+
+        tileRef['t_cols_offset'] = int(self.tile_cols * self.i_x)
+        tileRef['t_rows_offset'] = int(self.tile_rows * self.i_y)
+
+        if self.i_y >= self.num_tiles_y - 1:
+            self.i_y = 0
+            self.i_x += 1
+            # set bbox for next tile (E, W)      
+            self.tile_bbox['maxx'] += self.tile_length_x 
+            self.tile_bbox['minx'] += self.tile_length_x 
+        else:
+            self.i_y += 1
+
+        return query_url, tileRef
+
+    def _flipBbox(self, bbox, proj, version):
+        """ 
+        Flips coordinates if WMS standard is 1.3.0 and 
+        projection is geographic.
+
+        value flips between this keys:
+        maxy -> maxx
+        maxx -> maxy
+        miny -> minx
+        minx -> miny
+        @return copy of bbox with flipped coordinates
+        """  
+
+        temp_bbox = dict(bbox)
+        new_bbox = {}
+        new_bbox['maxy'] = temp_bbox['maxx']
+        new_bbox['miny'] = temp_bbox['minx']
+        new_bbox['maxx'] = temp_bbox['maxy']
+        new_bbox['minx'] = temp_bbox['miny']
+
+        return new_bbox
+
+
+class WMTSRequestMgr(BaseRequestMgr):
+    def __init__(self, params, bbox, region, proj_srs, cap_file = None):
+        """!Initializes data needed for iteration through tiles.
+        """
+
+        self.proj_srs = proj_srs
+        self.meters_per_unit = None
+
+        # constant defined in WMTS standard (in meters)
+        self.pixel_size = 0.00028
+
+        # parse capabilities file
+        cap_tree = etree.parse(cap_file)
+        root = cap_tree.getroot()
+
+        # get layer tile matrix sets with same projection
+        mat_sets = self._getMatSets(root, params)
+
+        # TODO: what if more tile matrix sets are returned?
+        params['tile_matrix_set'] = mat_sets[0].find(self._ns_ows('Identifier')).text
+
+        # find tile matrix with resolution closest and smaller to wanted resolution 
+        tile_mat  = self._findTileMats(mat_sets[0].findall(self._ns_wmts('TileMatrix')), region, bbox)
+
+        # initialize data needed for iteration through tiles
+        self._computeRequestData(tile_mat, params, bbox)
+
+    def GetMapRegion(self):
+        """!Get size in pixels and bounding box of raster where all tiles will be merged.
+        """
+        return self.map_region
+
+    def _getMatSets(self, root, params):
+        """!Get matrix sets which are available for chosen layer and have required EPSG.
+        """
+        contents = self.find(root, 'Contents', self._ns_wmts)
+        layers = self.findall(contents, 'Layer', self._ns_wmts)
+
+        ch_layer = None
+        for layer in layers:
+            layer_id = layer.find(self._ns_ows('Identifier')).text
+            if layer_id and layer_id == params['layers']:
+                ch_layer = layer 
+                break  
+
+        if ch_layer is None:
+            grass.fatal(_("Layer '%s' was not found in capabilities file") % params['layers'])
+
+        mat_set_links = self.findall(ch_layer, 'TileMatrixSetLink', self._ns_wmts)
+
+        suitable_mat_sets = []
+        tileMatrixSets = self.findall(contents, 'TileMatrixSet', self._ns_wmts)
+
+        for link in  mat_set_links:
+            mat_set_link_id = link.find(self._ns_wmts('TileMatrixSet')).text
+            if not mat_set_link_id:
+                continue
+            for mat_set in tileMatrixSets:
+                mat_set_id = mat_set.find(self._ns_ows('Identifier')).text 
+                if mat_set_id and mat_set_id != mat_set_link_id:
+                    continue
+                mat_set_srs = mat_set.find(self._ns_ows('SupportedCRS')).text
+                if mat_set_srs and mat_set_srs.lower() == ("EPSG:"+ str(params['srs'])).lower():
+                    suitable_mat_sets.append(mat_set)
+
+        if not suitable_mat_sets:
+            grass.fatal(_("Layer '%s' is not available with %s code.") % (params['layers'],  "EPSG:" + str(params['srs'])))
+
+        return suitable_mat_sets
+
+    def _findTileMats(self, tile_mats, region, bbox):
+        """!Find tile matrix set with closest and smaller resolution to requested resolution
+        """        
+        scale_dens = []
+
+        scale_dens.append((bbox['maxy'] - bbox['miny']) / region['rows'] * self._getMetersPerUnit()  / self.pixel_size)
+        scale_dens.append((bbox['maxx'] - bbox['minx']) / region['cols'] * self._getMetersPerUnit() / self.pixel_size)
+
+        scale_den = min(scale_dens)
+
+        first = True
+        for t_mat in tile_mats:
+            mat_scale_den = float(t_mat.find(self._ns_wmts('ScaleDenominator')).text)
+            if first:
+                best_scale_den = mat_scale_den
+                best_t_mat = t_mat
+                first = False
+                continue
+                
+            best_diff = best_scale_den - scale_den
+            mat_diff = mat_scale_den - scale_den
+            if (best_diff < mat_diff  and  mat_diff < 0) or \
+               (best_diff > mat_diff  and  best_diff > 0):
+                best_t_mat = t_mat
+                best_scale_den = mat_scale_den
+
+        return best_t_mat
+
+    def _getMetersPerUnit(self):
+        """!Get coefficient which allows to convert units of request projection into meters 
+        """  
+        if self.meters_per_unit:
+            return self.meters_per_unit
+
+        # for geographic projection
+        if self._isGeoProj(self.proj_srs):
+            proj_params = proj_srs.split(' ')
+            for param in proj_parmas:
+                if '+a' in param:
+                    a = float(param.split('=')[1])
+                    break
+            equator_perim = 2 * pi * a
+            # meters per degree on equator
+            self.meters_per_unit = equator_perim / 360 
+
+        # other units
+        elif '+to_meter' in self.proj_srs:
+            proj_params = self.proj_srs.split(' ')
+            for param in proj_params:
+                if '+to_meter' in param:
+                    self.meters_per_unit = 1/float(param.split('=')[1])
+                    break
+        # coordinate system in meters        
+        else:
+            self.meters_per_unit = 1
+
+        return self.meters_per_unit
+
+    def _computeRequestData(self, tile_mat, params, bbox):
+        """!Initialize data needed for iteration through tiles.
+        """  
+        scale_den = float(tile_mat.find(self._ns_wmts('ScaleDenominator')).text)
+
+        pixel_span = scale_den * self.pixel_size / self._getMetersPerUnit()
+
+        tl_corner = tile_mat.find(self._ns_wmts('TopLeftCorner')).text.split(' ')
+
+        self.t_cols = int(tile_mat.find(self._ns_wmts('TileWidth')).text)
+        t_length_x = pixel_span * self.t_cols
+        
+        self.t_rows = int(tile_mat.find(self._ns_wmts('TileHeight')).text)
+        t_length_y = pixel_span * self.t_rows
+    
+        t_mat_min_x = float(tl_corner[0])
+        t_mat_max_y = float(tl_corner[1])
+        
+        epsilon = 2e-15
+
+        self.t_num_bbox = {}
+        self.t_num_bbox['min_col'] = int(floor((bbox['minx'] - t_mat_min_x) / t_length_x + epsilon))
+        self.t_num_bbox['max_col'] = int(floor((bbox['maxx'] - t_mat_min_x) / t_length_x - epsilon))
+
+        self.t_num_bbox['min_row'] = int(floor((t_mat_max_y - bbox['maxy']) / t_length_y + epsilon))
+        self.t_num_bbox['max_row'] = int(floor((t_mat_max_y - bbox['miny']) / t_length_y - epsilon))
+
+        matWidth = int(tile_mat.find(self._ns_wmts('MatrixWidth')).text)
+        matHeight = int(tile_mat.find(self._ns_wmts('MatrixHeight')).text)
+
+        self.intersects = False
+        for col in ['min_col', 'max_col']:
+            for row in ['min_row', 'max_row']:
+                if (0 <= self.t_num_bbox[row] and self.t_num_bbox[row] <= matHeight) and \
+                   (0 <= self.t_num_bbox[col] and self.t_num_bbox[col] <= matWidth):
+                    self.intersects = True
+                    
+        if not self.intersects:
+            grass.warning(_('Region is out of server data extend.'))
+            self.map_region = None
+            return
+
+        if self.t_num_bbox['min_col'] < 0:
+            self.t_num_bbox['min_col']  = 0
+
+        if self.t_num_bbox['max_col'] > (matWidth - 1):
+            self.t_num_bbox['max_col']  = matWidth - 1
+
+        if self.t_num_bbox['min_row'] < 0:
+            self.t_num_bbox['min_row']  = 0
+
+        if self.t_num_bbox['max_row'] > (matHeight - 1):
+            self.t_num_bbox['max_row'] = matHeight - 1
+
+        num_tiles = (self.t_num_bbox['max_col'] - self.t_num_bbox['min_col'] + 1) * (self.t_num_bbox['max_row'] - self.t_num_bbox['min_row'] + 1) 
+        grass.message(_('Fetching %d tiles with %d x %d pixel size per tile...') % (num_tiles, self.t_cols, self.t_rows))
+
+        self.url = params['url'] + ("?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&" \
+                                    "LAYER=%s&STYLE=%s&FORMAT=%s&TILEMATRIXSET=%s&TILEMATRIX=%s" % \
+                                   (params['layers'], params['styles'], params['format'],
+                                    params['tile_matrix_set'], tile_mat.find(self._ns_ows('Identifier')).text ))
+
+        if params['urlparams'] != "":
+            self.url += "&" + params['urlparams']
+
+        self.map_region = {}
+        self.map_region['minx'] = self.t_num_bbox['min_col'] * t_length_x + t_mat_min_x
+        self.map_region['maxy'] = t_mat_max_y - (self.t_num_bbox['min_row']) * t_length_y
+
+        self.map_region['maxx'] = (self.t_num_bbox['max_col'] + 1) * t_length_x + t_mat_min_x
+        self.map_region['miny'] = t_mat_max_y - (self.t_num_bbox['max_row'] + 1) * t_length_y
+
+        self.map_region['cols'] = self.t_cols * (self.t_num_bbox['max_col'] -  self.t_num_bbox['min_col'] + 1)
+        self.map_region['rows'] = self.t_rows * (self.t_num_bbox['max_row'] -  self.t_num_bbox['min_row'] + 1)
+
+        self.i_col = self.t_num_bbox['min_col']
+        self.i_row = self.t_num_bbox['min_row'] 
+
+        self.tileRef = {
+                          'sizeX' : self.t_cols,
+                          'sizeY' : self.t_rows
+                        }
+
+    def GetNextTile(self):
+        """!Get url for request the tile from server and information for merging the tile with other tiles
+        """
+        if not self.intersects or self.i_col > self.t_num_bbox['max_col']:
+            return None 
+                
+        query_url = self.url + "&TILECOL=%i&TILEROW=%i" % (int(self.i_col), int(self.i_row)) 
+
+        self.tileRef['t_cols_offset'] = int(self.t_cols * (self.i_col - self.t_num_bbox['min_col']))
+        self.tileRef['t_rows_offset'] = int(self.t_rows * (self.i_row - self.t_num_bbox['min_row']))
+
+        if self.i_row >= self.t_num_bbox['max_row']:
+            self.i_row =  self.t_num_bbox['min_row']
+            self.i_col += 1
+        else:
+            self.i_row += 1 
+
+        return query_url, self.tileRef
+
+    def _ns_wmts(self, tag):
+        """!Helper method - XML namespace
+        """       
+        return "{http://www.opengis.net/wmts/1.0}" + tag
+
+    def _ns_ows(self, tag):
+        """!Helper method - XML namespace
+        """
+        return "{http://www.opengis.net/ows/1.1}" + tag

Modified: grass-addons/grass7/raster/r.in.wms2/wms_gdal_drv.py
===================================================================
--- grass-addons/grass7/raster/r.in.wms2/wms_gdal_drv.py	2012-10-09 21:50:21 UTC (rev 53347)
+++ grass-addons/grass7/raster/r.in.wms2/wms_gdal_drv.py	2012-10-09 22:02:55 UTC (rev 53348)
@@ -1,4 +1,5 @@
 import os
+import grass.script as grass 
 
 try:
     from osgeo import gdal
@@ -8,8 +9,6 @@
 
 import xml.etree.ElementTree as etree
 
-import grass.script as grass
-
 from wms_base import WMSBase
 
 class NullDevice():
@@ -30,25 +29,25 @@
         service.set("name","WMS")
         
         version = etree.SubElement(service, "Version")
-        version.text =self.o_wms_version
+        version.text =self.params['wms_version']
         
         server_url = etree.SubElement(service, "ServerUrl")
-        server_url.text =self.o_mapserver_url
+        server_url.text =self.params['url']
         
-        srs = etree.SubElement(service, self.projection_name)   
-        srs.text = 'EPSG:' + str(self.o_srs)
+        srs = etree.SubElement(service, self.params['proj_name'])   
+        srs.text = 'EPSG:' + str(self.params['srs'])
         
         image_format = etree.SubElement(service, "ImageFormat")
-        image_format.text = self.mime_format
+        image_format.text = self.params['format']
         
         image_format = etree.SubElement(service, "Transparent")
-        image_format.text = self.transparent
+        image_format.text = self.params['transparent']
         
         layers = etree.SubElement(service, "Layers")
-        layers.text = self.o_layers
+        layers.text = self.params['layers']
         
         styles = etree.SubElement(service, "Styles")
-        styles.text = self.o_styles
+        styles.text = self.params['styles']
         
         data_window = etree.SubElement(gdal_wms, "DataWindow")
         
@@ -76,10 +75,10 @@
         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_x.text = str(self.tile_size['cols']) 
         
         block_size_y = etree.SubElement(gdal_wms, "BlockSizeY")
-        block_size_y.text = str(self.tile_rows)
+        block_size_y.text = str(self.tile_size['rows'])
         
         xml_file = self._tempfile()
         
@@ -98,9 +97,11 @@
 
         # GDAL WMS driver does not flip geographic coordinates 
         # according to WMS standard 1.3.0.
-        if self.flip_coords and self.o_wms_version == "1.3.0":
+        if ("+proj=latlong" in self.proj_srs or \
+            "+proj=longlat" in self.proj_srs) and \
+            self.params['wms_version'] == "1.3.0":
             grass.warning(_("If module will not be able to fetch the data in this\
-                           geographic projection, \n try flag -d or use WMS version 1.1.1."))
+                           geographic projection, \n try 'WMS_GRASS' driver or use WMS version 1.1.1."))
 
         self._debug("_download", "started")
         



More information about the grass-commit mailing list