Changes to allow for less strict testing of extent, 2D NetCDF rasters, smart...

Changes to allow for less strict testing of extent, 2D NetCDF rasters, smart searching of a unique value in a raster etc.

formats/asciigrid.py

@@ -33,7 +33,7 @@ class AsciiGrid(Raster, GridEnvelope2D):
         # Module wide constants
         self._const = Const()
         self._const.FILEXT = "asc";     
-        self._const.MAXDIGITSPERCELL = 6 # TODO this is hardcoded - change this
+        self._const.MAXDIGITSPERCELL = 8 # TODO this is hardcoded - change this
         self.name = "dummy." + self._const.FILEXT;
         
         # Initialise further

formats/gridenvelope2d.py

@@ -84,20 +84,25 @@ class GridEnvelope2D(object):
             cy = self.getMaxY() - (i + 0.5) * self.dy;
         return (cx, cy);
     
-    def hasSameExtent(self, obj):
+    def hasSameExtent(self, obj, epsilon=const.epsilon):
+        # Check input
         if (not isinstance(obj, GridEnvelope2D)): 
             return False;
+        if epsilon < 0:
+            return False;
+        
+        # Now check the extent
         if self.ncols != obj.ncols:
             return False;
         if self.nrows != obj.nrows:
             return False;
-        if abs(self.dy - obj.dy) > const.epsilon:
+        if abs(self.dy - obj.dy) > epsilon:
             return False;
-        if abs(self.dx - obj.dx) > const.epsilon:
+        if abs(self.dx - obj.dx) > epsilon:
             return False;
-        if abs(self.xll - obj.xll) > const.epsilon:
+        if abs(self.xll - obj.xll) > epsilon:
             return False;
-        if abs(self.yll - obj.yll) > const.epsilon:
+        if abs(self.yll - obj.yll) > epsilon:
             return False;
         else:
             return True;        

formats/netcdf4envelope2d.py

@@ -2,16 +2,16 @@ from .gridenvelope2d import GridEnvelope2D;
 
 class Netcdf4Envelope2D(GridEnvelope2D):
     # Constants
-    LON = 'lon';
-    LAT = 'lat';
+    __X = 'lon';
+    __Y = 'lat';
     TIME = 'time'
     
     # TODO: here the names for the dimensions are hardcoded, but they can be retrieved
     # from the netCDF4 file 
-    def __init__(self, ds):
+    def __init__(self, ds, X='lon', Y='lat'):
         # Initialise
-        LON = self.LON;
-        LAT = self.LAT;
+        self.__X = X;
+        self.__Y = Y;
         
         # Make sure you can call the __init__ method of the base class
         if ds != None:
@@ -19,19 +19,19 @@ class Netcdf4Envelope2D(GridEnvelope2D):
             self.__dataset = ds;
             _dims = ds.dimensions;
             _vars = ds.variables;
-            x_range= self._readRange(_vars[LON]);
-            y_range = self._readRange(_vars[LAT]);
+            x_range = self._readRange(_vars[self.__X]);
+            y_range = self._readRange(_vars[self.__Y]);
             
             # Now get the necessary properties for initialising an envelope object
-            dx = GridEnvelope2D._getStep(x_range[0], x_range[1], len(_dims[LON]));
-            dy = GridEnvelope2D._getStep(y_range[0], y_range[1], len(_dims[LAT]));
-            xll = min(_vars[LON]) - 0.5*dx;
-            yll = min(_vars[LAT]) - 0.5*dy;
-            GridEnvelope2D.__init__(self, len(_dims[LON]), len(_dims[LAT]), xll, yll, dx, dy);
+            dx = GridEnvelope2D._getStep(x_range[0], x_range[1], len(_dims[self.__X]));
+            dy = GridEnvelope2D._getStep(y_range[0], y_range[1], len(_dims[self.__Y]));
+            xll = min(_vars[self.__X]) - 0.5*dx;
+            yll = min(_vars[self.__Y]) - 0.5*dy;
+            GridEnvelope2D.__init__(self, len(_dims[self.__X]), len(_dims[self.__Y]), xll, yll, dx, dy);
         
             # In some netCDF files the longitudes and latitudes are stored differently
-            if _vars[LON][0] > _vars[LON][-1]: self.xcoords_sort = 'DESC';
-            if _vars[LAT][0] < _vars[LAT][-1]: self.ycoords_sort = 'ASC';
+            if _vars[self.__X][0] > _vars[self.__X][-1]: self.xcoords_sort = 'DESC';
+            if _vars[self.__Y][0] < _vars[self.__Y][-1]: self.ycoords_sort = 'ASC';
     
     def _readRange(self, varxy):
         # Argument is either x or y

formats/netcdf4raster.py

@@ -6,7 +6,7 @@ import os;
 class Netcdf4Raster(Netcdf4Envelope2D):
     # Constants
     DATAFILEXT = 'nc4';
-    _original_name = "yield_mai"
+    _original_name = "Band1" # TODO: set it when a template file is opened for appending
     
     # Data attributes - assign some dummy values for the mean time
     name = "dummy.nc4";
@@ -16,12 +16,18 @@ class Netcdf4Raster(Netcdf4Envelope2D):
     _currow = 0;
     cellsize = 1;
     nodatavalue = -9999.0;
+    __X = ""
+    __Y = ""
     
     def __init__(self, filepath):
         # Retrieve the name from the filepath and assign - incl. extension
         self.name = os.path.basename(filepath);
         # Also derive the folder
         self.folder = os.path.dirname(filepath);
+        
+    def set_xy_ranges(self, xvar, yvar):
+        self.__X = xvar
+        self.__Y = yvar
 
     def open(self, mode, ncols=1, nrows=1, xll=0, yll=0, cellsize=1, nodatavalue=-9999.0):
         # If file does not exist and mode[0] = 'w', create it!
@@ -30,13 +36,20 @@ class Netcdf4Raster(Netcdf4Envelope2D):
             if os.path.exists(fpath):
                 print("About to open file " + fpath + " in append mode");
             self._dataset = Dataset(fpath, 'a', format='NETCDF4');
-            Netcdf4Envelope2D.__init__(self, self._dataset);
+            self.nodatavalue = nodatavalue
+            if (self.__X == "") or (self.__Y == ""):
+                Netcdf4Envelope2D.__init__(self, self._dataset);
+            else:
+                Netcdf4Envelope2D.__init__(self, self._dataset, self.__X, self.__Y);
             return True;
         else:
             if os.path.exists(fpath):  
                 # Open the netCDF4 file          
                 ds = Dataset(fpath, 'r');
-                Netcdf4Envelope2D.__init__(self, ds);
+                if (self.__X == "") or (self.__Y == ""):
+                    Netcdf4Envelope2D.__init__(self, ds);
+                else:
+                    Netcdf4Envelope2D.__init__(self, self._dataset, self.__X, self.__Y);   
                 self._dataset = ds;
                 self._varname = self._get_varname();
                 self.nodatavalue = ds.variables[self._varname]._FillValue
@@ -65,7 +78,8 @@ class Netcdf4Raster(Netcdf4Envelope2D):
         if self._varname != "":
             if parseLine:
                 # Indexing: 1. years 2. by rows 3. columns
-                result = self.getVariables(self._varname)[:, self._currow, :]; 
+                # 3D: result = self.getVariables(self._varname)[:, self._currow, :]; 
+                result = self.getVariables(self._varname)[self._currow, :] # 2D
             self._currow += 1;  # row index is zero-based! 
             return result;
     
@@ -74,12 +88,14 @@ class Netcdf4Raster(Netcdf4Envelope2D):
         return self.DATAFILEXT;
     
     def writeheader(self, name, long_name, units):
+        # TODO: if _original_name is not one of the variables in the file -> error
         v = self._dataset.variables[self._original_name]
         v.setncattr("long_name", long_name)
         v.setncattr("units", units)
         if self._original_name != name:
             self._dataset.renameVariable(self._original_name, name)
         self._varname = name; 
+        self._dataset.variables[self._varname].missing_value = self.nodatavalue
     
     def writenext(self, sequence_with_data):
         # Assume that the sequence is indexed 1. by year and 2. by column
@@ -90,12 +106,15 @@ class Netcdf4Raster(Netcdf4Envelope2D):
             raise ValueError("Input array has unexpected shape")
         if sequence_with_data.shape[1] != self.ncols:
             raise ValueError("Input array has unexpected dimension")
-        self._dataset.variables[key][:, self._currow, :] = sequence_with_data
+        # TODO adapt so that Netcdf files with different dimensions can be written
+        # 3D: self._dataset.variables[key][:, self._currow, :] = sequence_with_data
+        self._dataset.variables[key][self._currow, :] = sequence_with_data #2D
         self._currow += 1;
         
     def close(self):
-        if self._dataset:
-            self._dataset.close(); 
+        if not self._dataset is None:
+            if self._dataset.isopen():
+                self._dataset.close(); 
 
     def reset(self):
         self._currow = 0; 
@@ -114,6 +133,13 @@ class Netcdf4Raster(Netcdf4Envelope2D):
     def getVariableName(self):
         return self._varname;
     
+    def get_value(self, i, k):    
+        # Return the wanted value
+        for _ in range(0, i): self.next(False)
+        line = self.next()
+        self.reset()
+        return line[int(k)]
+    
 """
 def write_nc4_file(filename, crop_no):
     # Write the data to the file

toolbox/cliplib.py

 """ClipLib - a Python library for clipping from rasters"""
-from formats.gridenvelope2d import GridEnvelope2D
+from lmgeo.formats.gridenvelope2d import GridEnvelope2D
 from array import array
 try:
     import numpy as np

toolbox/punchlib.py

 import numpy as np
-import toolbox.cliplib as cp
-import formats.gridenvelope2d as ge
+import lmgeo.toolbox.cliplib as cp
+import lmgeo.formats.gridenvelope2d as ge
 from pyproj import Proj, transform
 
 # Constants
@@ -264,9 +264,37 @@ def get_zonal_stats_as_table(shapeRecords, srchfld, rastergrid, stat_types):
         result.append(rec)
     return result
 
-def get_centroid(rastergrid, i, k):
-    result = [rastergrid.xll + (k+0.5)*rastergrid.dx]
-    result.extend([rastergrid.yll + (rastergrid.nrows-i-0.5)*rastergrid.dy])
+def get_centroid(rastergrid, colidx, rowidx):
+    result = [rastergrid.xll + (colidx+0.5)*rastergrid.dx]
+    result.extend([rastergrid.yll + (rastergrid.nrows-rowidx-0.5)*rastergrid.dy])
     return result
-        
+
+def get_coords_from_id(rastergrid, id, startrow=0, endrow=-1):
+    # Initialise
+    result = None
+    if endrow == -1: endrow = rastergrid.nrows - 1
+    
+    # The rastergrid is assumed to hold unique values
+    for i in range(startrow): rastergrid.next(False)
+    for i in range(startrow, endrow+1):
+        line = rastergrid.next(True)
+        for k in range(rastergrid.ncols):
+            if line[k] == id:
+                result = get_centroid(rastergrid, k, i)
+                rastergrid.reset()
+                return result
+            else:
+
+                if (i == endrow) and (k == rastergrid.ncols - 1):
+                    if (startrow == 0):
+                        # We searched everywhere but did not find
+                        rastergrid.reset()
+                        return result # None
+                    else:
+                        # We have reached the end of the file w/o success
+                        rastergrid.reset()
+                        result = get_coords_from_id(rastergrid, id, endrow=startrow-1)
+                        return result
+                        
+    
     
\ No newline at end of file