Merge branch 'develop'

.gitignore

@@ -7,4 +7,9 @@ MANIFEST
 *.tar.*
 build
 dist
-hex_utils/__pycache__
+hex_utils/__pycache__
+movies
+data
+docs
+tests
+helpers

.pydevproject

@@ -2,6 +2,7 @@
 <?eclipse-pydev version="1.0"?><pydev_project>
 <pydev_pathproperty name="org.python.pydev.PROJECT_SOURCE_PATH">
 <path>/${PROJECT_DIR_NAME}</path>
+<path>/${PROJECT_DIR_NAME}/helpers</path>
 </pydev_pathproperty>
 <pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 3.0</pydev_property>
 <pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Default</pydev_property>

hex_utils/asc.py

@@ -10,7 +10,8 @@
 #
 # [0] http://resources.esri.com/help/9.3/arcgisengine/java/GP_ToolRef/spatial_analyst_tools/esri_ascii_raster_format.htm
 
-import math
+import sys, math
+from scipy import interpolate
 from hex_utils.grid import Grid
 
 class ASC (Grid):
@@ -42,7 +43,7 @@ def init(self, ncols, nrows, xll, yll, size, nodata = ""):
         self._set_size(size)
 
 
-    def getNearestNeighbour(self, x, y):
+    def _getNearestNeighbourGridCoords(self, x, y):
 
         if x < self._xll:
             x = self._xll
@@ -58,14 +59,84 @@ def getNearestNeighbour(self, x, y):
 
         i = math.trunc((x - self._xll) / self._size)
         j = self._nrows - 1 - math.trunc((y - self._yll) / self._size)
+
+        return i, j
+
+
+    def getNearestNeighbour(self, x, y):
+
+        i, j = self._getNearestNeighbourGridCoords(x, y)
 
         try:
             return self._grid[i][j]
         except IndexError:
             raise IndexError("Wrong indexes in nearest neighbour:" + 
                 "i: " + str(i) + " j: " + str(j) + " x: " + str(x) + " y: " + str(y))
+
+
+    def _getNeighbourhoodGridCoords(self, i, j):
+
+        ii = []
+        jj = []
+
+        if i > 0:
+
+            if j > 0:
+                ii.append(i-1)
+                jj.append(j-1)
+
+            ii.append(i-1)
+            jj.append(j)
+
+            if j < self.nrows - 1:
+                ii.append(i-1)
+                jj.append(j+1)
+
+        if j > 0:
+            ii.append(i)
+            jj.append(j-1)
+
+        ii.append(i)
+        jj.append(j)
+
+        if j < self.nrows - 1:
+            ii.append(i)
+            jj.append(j+1)
+
+        if i < self.ncols - 1:
+
+            if j > 0:
+                ii.append(i+1)
+                jj.append(j-1)
+
+            ii.append(i+1)
+            jj.append(j)
+
+            if j < self.nrows - 1:
+                ii.append(i+1)
+                jj.append(j+1)
+
+        return ii, jj
+
 
+    def interpolMultiquadratic(self, x, y, epsilon=100):
+
+        xx = []
+        yy = []
+        vals = []
+        i, j = self._getNearestNeighbourGridCoords(x, y)
+        ii, jj = self._getNeighbourhoodGridCoords(i, j)
+
+        for n in range(len(ii)):
+            if ii[n] != None and jj[n] != None and self._grid[ii[n]][jj[n]] != None:
+                xx.append(self._xll + ii[n] * self._size + self._size / 2)
+                yy.append(self._yll + self._nrows * self._size - jj[n] * self._size - self._size / 2)
+                vals.append(self._grid[ii[n]][jj[n]])
+
+        f = interpolate.Rbf(xx, yy, vals, epsilon=epsilon)
+        return f(x,y)
 
+
     def _loadHeader(self):
 
         # Mandatory header

hex_utils/asc2hasc.py

@@ -6,8 +6,8 @@
 # Transforms rectangular ESRI ASCCI grids into ASCII encoded cartographic 
 # hexagonal grids (HASC) [0]. The resulting hexagonal grid is created with a
 # spatial resolution similar to the input rectangular grid. The values of the
-# resulting hexagonal grid are computed using a simple nearest-neighbour 
-# algorithm.
+# resulting hexagonal grid are computed using a simple nearest-neighbour or
+# a multi-quadratic interpolation.
 #
 # Author: Luís Moreira de Sousa (luis.de.sousa[@]protonmail.ch)
 # Date: 31-03-2016 
@@ -18,52 +18,55 @@
 import math
 from hex_utils.asc import ASC
 from hex_utils.hasc import HASC
+import argparse 
+from enum import Enum  
 
-RES_FACTOR = 1.134
+RES_FACTOR = 1.134   
 
-def wrongUsage():
-
-    print("\nThis programme requires three arguments:\n" +
-          " - conversion mode \n" +
-          " - path to an input ESRI ASCII file \n" +
-          " - path to the output HASC file \n\n" + 
-          "The conversion can be of two types: \n" +
-          " -a : preserve cell area  \n" +
-          " -r : preserve spatial resolution \n\n" +
-          "Usage example: \n"
-          "   asc2hasc -r /path/to/input.asc /path/to/output.hasc\n")
-    sys.exit()
+class Method(Enum):
+    MULTIQUADRATIC = 'mq',
+    NEAREST_NEIGHBOUR = 'nn' 
 
-# This method is the same as in other command line utils
-def processArguments(args):
+
+
+def setArguments():
 
-    global inputFile
-    global outputFile
+    parser = argparse.ArgumentParser(description='Converts an ESRI ASCII grid into an HexASCII (HASC) grid.')
+    parser.add_argument("-a", "--area", dest="area", default = 0,
+                      type=float, help="cell area of the output grid" )
+    parser.add_argument("-r", "--resolution", action='store_true',
+                      help = "preserve original spatial resolution")
+    parser.add_argument("-m", "--method", default="mq",
+                      help = "Interpolation method: mq - Multiquadratic, nn - Nearest Neighbour")
+    parser.add_argument("-i", "--input", dest="inputFile", required = True,
+                      help="input ESRI ASCII grid file" )
+    parser.add_argument("-o", "--output", dest="outputFile", default = "out.hasc",
+                      help="output HexASCII grid file" )
 
-    if len(args) < 4 or (str(args[1]) != '-r' and str(args[1]) != '-a'):
-        wrongUsage() 
-    else:
-        inputFile = str(args[2])
-        outputFile = str(args[3])
+    return parser.parse_args()
+
 
 # ------------ Main ------------ #
 def main():
 
-    processArguments(sys.argv)
-     
+    args = setArguments()
+
     esriGrid = ASC()
     try:
-        esriGrid.loadFromFile(inputFile)
+        esriGrid.loadFromFile(args.inputFile)
     except (ValueError, IOError) as ex:
-        print("Error loading the grid %s: %s" % (inputFile, ex))
+        print("Error loading the grid %s: %s" % (args.inputFile, ex))
         sys.exit()
 
     esriArea = math.pow(esriGrid.size, 2)
     hexArea = 0
 
     # Preserve spatial resolution: increase cell area.
-    if (sys.argv[1] == '-r'): 
+    if (args.resolution): 
         hexArea = esriArea * RES_FACTOR
+    # Use area provided.
+    elif (args.area > 0):
+        hexArea = args.area
     # Preserve cell area: used the same as the original grid.
     else:
         hexArea = esriArea
@@ -96,12 +99,17 @@ def main():
     hexGrid = HASC()
     hexGrid.init(hexCols, hexRows, hexXLL, hexYLL, hexSide, esriGrid.nodata)
 
+    if(args.method == Method.MULTIQUADRATIC):
+        interpol = esriGrid.getNearestNeighbour
+    else:
+        interpol = esriGrid.interpolMultiquadratic    
+
     for j in range(hexGrid.nrows):
         for i in range(hexGrid.ncols):
             x, y = hexGrid.getCellCentroidCoords(i, j)
-            hexGrid.set(i, j, esriGrid.getNearestNeighbour(x, y))
+            hexGrid.set(i, j, interpol(x, y))
 
-    hexGrid.save(outputFile)
+    hexGrid.save(args.outputFile)
 
     print ("Finished successfully.")
 

hex_utils/grid.py

@@ -148,10 +148,10 @@ def _loadLineValues(self, values):
                     "Expected something in the range [0.." + 
                     str(self._ncols - 1) + ", 0.." + str(self._nrows - 1) + "]")
 
-            self._colIdx += 1;
+            self._colIdx += 1
             if self._colIdx >= self._ncols:
-                self._colIdx = 0;
-                self._rowIdx += 1;               
+                self._colIdx = 0
+                self._rowIdx += 1               
 
 
     def _loadValues(self):
@@ -167,7 +167,7 @@ def _loadValues(self):
         while (self._nextLine):
             self._loadLineValues(self._nextLine.split())
             self._nextLine = self._file.readline()
-
+            
 
     def _saveHeader(self, f):
         raise NotImplementedError("Please Implement this method")

requirements.txt

@@ -0,0 +1,11 @@
+cycler==0.10.0
+GDAL==2.1.0
+hex-utils==0.2
+matplotlib==1.5.3
+numpy==1.11.1
+pkg-resources==0.0.0
+pyparsing==2.1.9
+python-dateutil==2.5.3
+pytz==2016.6.1
+scipy==0.18.1
+six==1.10.0

setup.py

@@ -12,8 +12,8 @@
         'console_scripts': [
             'hasc2gml=hex_utils.hasc2gml:main',
             'asc2hasc=hex_utils.asc2hasc:main',
-            'asc2hasc=hex_utils.surface2hasc:main',
-            'asc2hasc=hex_utils.surface2asc:main'
+            'surface2hasc=hex_utils.surface2hasc:main',
+            'surface2asc=hex_utils.surface2asc:main'
         ],
     },
 
@@ -36,4 +36,4 @@
         ],
 
     # could also include long_description, download_url, classifiers, etc.
-)
+)

surfaces/surfaceGaussian.py

@@ -37,5 +37,5 @@ def fun(x, y):
     return gauss.fun(x, y)
 
 # Uncomment this lines for auto-plotting
-gauss.plot()
+# gauss.plot()
 

surfaces/surfaceHubbert.py

@@ -39,4 +39,4 @@ def fun(x, y):
     return hubbert.fun(x, y)
 
 # Uncomment these lines for auto-plotting
-hubbert.plot()
+# hubbert.plot()