Made resolution one number (isotropic) instead of two and added it to…

… keywords when getting layer metadata. This is in preparation for issues #168. Came across new issue enroute: #170

impact/storage/io.py

@@ -230,19 +230,23 @@ def get_metadata_from_layer(layer):
     if layer.datatype == 'raster':
         geotransform = extract_geotransform(layer)
         metadata['geotransform'] = geotransform
-        metadata['resolution'] = geotransform2resolution(geotransform)
+        metadata['resolution'] = geotransform2resolution(geotransform,
+                                                         isotropic=True)
+    else:
+        metadata['resolution'] = None
+        metadata['geotransform'] = None
 
     # Metadata common to both raster and vector data
     metadata['bounding_box'] = layer.boundingBoxWGS84
     metadata['title'] = layer.title  # This maybe overwritten by keyword
     metadata['id'] = layer.id
 
     # Extract keywords
+    keyword_dict = {}
     if not hasattr(layer, 'keywords'):
-        msg = 'No keywords in %s. Submit patch to OWSLib maintainers' % layer
+        msg = 'No keywords in %s. Submit patch to OWSLib maintainers.' % layer
         raise Exception(msg)
     else:
-        keyword_dict = {}
         for keyword in layer.keywords:
             if keyword is not None:
                 # FIXME (Ole): Why would this be None sometimes?
@@ -254,11 +258,16 @@ def get_metadata_from_layer(layer):
                     else:
                         keyword_dict[keyword_string] = None
 
-        # Add a few extras:
-        # layertype
-        # resolution
-        metadata['keywords'] = keyword_dict
+    # Add resolution (as a string) so that layer "remembers"
+    # its original resolution
+    keyword_dict['resolution'] = str(metadata['resolution'])
+
+    # FIXME (Ole): This does not cause an Exception, and nothing is written to the log file
+    #              See issue #170
+    #raise Exception('weird')
 
+    # Record all keywords as part of the metadata and return
+    metadata['keywords'] = keyword_dict
     return metadata
 
 
@@ -289,7 +298,6 @@ def get_metadata(server_url, layer_name=None):
     # Get metadata for requested layer(s)
     metadata = {}
     for name in layer_names:
-
         if name in wcs.contents:
             layer = wcs.contents[name]
             layer.datatype = 'raster'  # Monkey patch layer type
@@ -541,6 +549,7 @@ def download(server_url, layer_name, bbox, resolution=None):
         if resolution is None:
             # Get native resolution and use that
             resolution = layer_metadata['resolution']
+            resolution = (resolution, resolution)  #FIXME (Ole): Make nicer
 
         # Download raster using specified bounding box and resolution
         template = WCS_TEMPLATE

impact/storage/raster.py

@@ -477,22 +477,20 @@ def get_resolution(self, isotropic=False):
         """Get raster resolution as a 2-tuple (resx, resy)
 
         Input
-            isotropic:
+            isotropic: If True, verify that dx == dy and return dx
+                       If False return 2-tuple (dx, dy)
         """
 
-        resx, resy = geotransform2resolution(self.geotransform)
-
-        if isotropic:
-            msg = ('Resolution for layer %s requested with '
-                   'isotropic=True, but '
-                   'resolutions in the horizontal and vertical '
-                   'are different: resx = %.12f, resy = %.12f. '
-                   % (self.get_name(), resx, resy))
-            assert numpy.allclose(resx, resy,
-                                  rtol=1.0e-6, atol=1.0e-8), msg
-            return resx
-        else:
-            return resx, resy
+        try:
+            res = geotransform2resolution(self.geotransform,
+                                          isotropic=isotropic)
+        except Exception, e:
+            msg = ('Resolution for layer %s could not be obtained: %s '
+                   % (self.get_name(), str(e)))
+            raise Exception(msg)
+
+        return res
+
 
     @property
     def is_raster(self):

impact/storage/utilities.py

@@ -173,16 +173,18 @@ def write_keywords(keywords, filename):
            'Expected %s.keywords' % (filename, basename))
     assert ext == '.keywords', msg
 
-
-
     # Write
     fid = open(filename, 'w')
     for k, v in keywords.items():
 
         msg = ('Key in keywords dictionary must be a string. '
-               'I got %s with type %s' % (k, type(k)))
+               'I got %s with type %s' % (k, str(type(k))[1:-1]))
         assert isinstance(k, basestring), msg
 
+        msg = ('Keyword value must be a string. '
+               'For key %s, I got %s with type %s' % (k, v, str(type(v))[1:-1]))
+        assert isinstance(v, basestring), msg
+
         key = k.strip()
 
         msg = ('Key in keywords dictionary must not contain the ":" '
@@ -321,25 +323,44 @@ def geotransform2bbox(geotransform, columns, rows):
     return [minx, miny, maxx, maxy]
 
 
-def geotransform2resolution(geotransform):
+def geotransform2resolution(geotransform, isotropic=False,
+                            # FIXME (Ole): Check these tolerances
+                            rtol=1.0e-2, atol=1.0e-5):
     """Convert geotransform to resolution
 
     Input
         geotransform: GDAL geotransform (6-tuple).
                       (top left x, w-e pixel resolution, rotation,
                       top left y, rotation, n-s pixel resolution).
                       See e.g. http://www.gdal.org/gdal_tutorial.html
+        Input
+            isotropic: If True, verify that dx == dy and return dx
+                       If False (default) return 2-tuple (dx, dy)
+            rtol, atol: Used to control how close dx and dy must be
+                        to quality for isotropic. These are passed on to
+                        numpy.allclose for comparison.
 
     Output
         resolution: grid spacing (resx, resy) in (positive) decimal
                     degrees ordered as longitude first, then latitude.
+                    or resx (if isotropic is True)
     """
 
-    x_res = geotransform[1]     # w-e pixel resolution
-    y_res = - geotransform[5]   # n-s pixel resolution (always negative)
+    resx = geotransform[1]     # w-e pixel resolution
+    resy = - geotransform[5]   # n-s pixel resolution (always negative)
 
-    return x_res, y_res
+    if isotropic:
+        msg = ('Resolution requested with '
+               'isotropic=True, but '
+               'resolutions in the horizontal and vertical '
+               'are different: resx = %.12f, resy = %.12f. '
+               % (resx, resy))
+        assert numpy.allclose(resx, resy,
+                              rtol=rtol, atol=atol), msg
 
+        return resx
+    else:
+        return resx, resy
 
 def bbox_intersection(*args):
     """Compute intersection between two or more bounding boxes
@@ -437,14 +458,14 @@ def buffered_bounding_box(bbox, resolution):
 
     Input
         bbox: Bounding box with format [W, S, E, N]
-        resolution: (resx, resy) Raster resolution in each direction.
+        resolution: (resx, resy) - Raster resolution in each direction.
+                    res - Raster resolution in either direction
                     If resolution is None bbox is returned unchanged.
 
     Ouput
         Adjusted bounding box
 
 
-
     Case in point: Interpolation point O would fall outside this domain
                    even though there are enough grid points to support it
 
@@ -462,20 +483,19 @@ def buffered_bounding_box(bbox, resolution):
     if resolution is None:
         return bbox
 
-    resx, resy = resolution
+    try:
+        resx, resy = resolution
+    except:
+        resx = resy = resolution
 
-    #print 'bbox', bbox
-    #print resx, resy
     bbox[0] -= resx
     bbox[1] -= resy
     bbox[2] += resx
     bbox[3] += resy
 
-    #print 'grown bbox', bbox
     return bbox
 
 
-
 def get_geometry_type(geometry):
     """Determine geometry type based on data
 

impact/tests/test_io.py

@@ -1403,12 +1403,15 @@ def test_geotransform2resolution(self):
         """
 
         for gt in GEOTRANSFORMS:
-            res = geotransform2resolution(gt)
-
+            res = geotransform2resolution(gt, isotropic=False)
             assert len(res) == 2
             assert numpy.allclose(res[0], gt[1], rtol=0, atol=1.0e-12)
             assert numpy.allclose(res[1], - gt[5], rtol=0, atol=1.0e-12)
 
+            res = geotransform2resolution(gt, isotropic=True)
+            assert numpy.allclose(res, gt[1], rtol=0, atol=1.0e-12)
+            assert numpy.allclose(res, - gt[5], rtol=0, atol=1.0e-12)
+
 if __name__ == '__main__':
     suite = unittest.makeSuite(Test_IO, 'test')
     runner = unittest.TextTestRunner(verbosity=2)

impact/views.py

@@ -126,10 +126,7 @@ def calculate(request, save_output=save_to_geonode):
             raster_resolution = None
         else:
             # Take the minimum
-            resx = min(haz_res[0], exp_res[0])
-            resy = min(haz_res[1], exp_res[1])
-
-            raster_resolution = (resx, resy)
+            raster_resolution = min(haz_res, exp_res)
 
         # New bounding box for data common to hazard, exposure and viewport
         # Download only data within this intersection
@@ -194,7 +191,12 @@ def calculate(request, save_output=save_to_geonode):
                              user=theuser)
     except Exception, e:
         # FIXME: Reimplement error saving for calculation.
-        # FIXME (Ole): Why?
+        # FIXME (Ole): Why should we reimplement?
+        # This is dangerous. Try to raise an exception
+        # e.g. in get_metadata_from_layer. Things will silently fail.
+        # See issue #170
+        #print 'ERRRORRRRRRR'
+
         logger.error(e)
         errors = e.__str__()
         trace = exception_format(e)