Merge branch 'dev'

geokit/__init__.py

@@ -40,7 +40,7 @@
 from geokit.core.raster import drawRaster
 
 # import the special algorithms
-#import geokit.algorithms
+import geokit.algorithms
 
 # Add useful paths for testing and stuff
 from collections import OrderedDict as _OrderedDict

geokit/core/extent.py

@@ -262,7 +262,7 @@ def __str__(s):
         return "(%.5f,%.5f,%.5f,%.5f)"%s.xyXY
 
 
-    def pad(s, pad): 
+    def pad(s, pad, percent=False): 
         """Pad the extent in all directions
         
         Parameters:
@@ -271,6 +271,10 @@ def pad(s, pad):
             The amount to pad in all directions
             * In units of the extent's srs
             * Can also accept a negative padding
+
+        percent : bool, optional
+            If True, the padding values are understood to be a percentage of the
+            unpadded extent
         
         Returns:
         --------
@@ -287,6 +291,10 @@ def pad(s, pad):
             xpad = pad
             ypad = pad
 
+        if percent:
+            xpad = xpad/100*(s.xMax-s.xMin) /2
+            ypad = ypad/100*(s.yMax-s.yMin) /2
+
         # Pad the extent
         return Extent(s.xMin-xpad, s.yMin-ypad, s.xMax+xpad, s.yMax+ypad, srs=s.srs)
 

geokit/core/location.py

@@ -258,6 +258,10 @@ def __init__(s, locations, srs=4326, _skip_check=False):
         if not _skip_check:
             if isinstance(locations, ogr.Geometry) or isinstance(locations, Location):
                 s._locations = np.array([Location.load(locations, srs=srs), ])
+            elif isinstance(locations, LocationSet):
+                s._locations = locations[:]
+            elif isinstance(locations, pd.DataFrame): 
+                s._locations = LocationSet( locations["geom"] )[:]
             else:
                 try: # Try loading all locations one at a time
                     s._locations = np.array([Location.load(l, srs=srs) for l in locations])
@@ -274,8 +278,11 @@ def __init__(s, locations, srs=4326, _skip_check=False):
         s._lats = None
         s._bounds4326 = None
         s.count = len(s._locations)
-
+        s.shape = (s.count,)
+
     def __len__(s): return s.count
+
+
     def __getitem__(s,i): return s._locations[i]
     def __repr__(s):
         out = " , Lon      , Lat\n"

geokit/core/raster.py

@@ -61,7 +61,7 @@ def gdalType(s):
 
 ####################################################################
 # Raster writer
-def createRaster( bounds, output=None, pixelWidth=100, pixelHeight=100, dtype=None, srs='europe_m', compress=True, noData=None, overwrite=True, fill=None, data=None, meta=None, **kwargs):
+def createRaster( bounds, output=None, pixelWidth=100, pixelHeight=100, dtype=None, srs='europe_m', compress=True, noData=None, overwrite=True, fill=None, data=None, meta=None, scale=1, offset=0, **kwargs):
     """Create a raster file
     
     NOTE:
@@ -135,6 +135,16 @@ def createRaster( bounds, output=None, pixelWidth=100, pixelHeight=100, dtype=No
         * array dimensions must fit raster dimensions as calculated by the bounds
           and the pixel resolution
 
+    scale : numeric; optional
+        The scaling value given to apply to all values
+        - numeric
+        * Must be the same datatye as the 'dtype' input (or that which is derived)
+
+    offset : numeric; optional
+        The offset value given to apply to all values
+        - numeric
+        * Must be the same datatye as the 'dtype' input (or that which is derived)
+
     Returns:
     --------
     * If 'output' is None: gdal.Dataset
@@ -152,7 +162,7 @@ def createRaster( bounds, output=None, pixelWidth=100, pixelHeight=100, dtype=No
                 raise GeoKitRasterError("Output file already exists: %s" %output)
 
     # Ensure bounds is okay
-    bounds = fitBoundsTo(bounds, pixelWidth, pixelHeight)
+    # bounds = fitBoundsTo(bounds, pixelWidth, pixelHeight)
 
     ## Make a raster dataset and pull the band/maskBand objects
     originX = bounds[0]
@@ -194,6 +204,8 @@ def createRaster( bounds, output=None, pixelWidth=100, pixelHeight=100, dtype=No
 
         # Fill the raster will zeros, null values, or initial values (if given)
         band = raster.GetRasterBand(1)
+        if not scale is None: band.SetScale(scale)
+        if not offset is None: band.SetOffset(offset)
 
         if( not noData is None):
             band.SetNoDataValue(noData)
@@ -367,8 +379,12 @@ def extractMatrix(source, bounds=None, boundsSRS='latlon', maskBand=False, autoc
         ywin=None
 
     # get Data
-    if maskBand: data = mb.ReadAsArray(xoff=xoff, yoff=yoff, win_xsize=xwin, win_ysize=ywin)
-    else: data = sourceBand.ReadAsArray(xoff=xoff, yoff=yoff, win_xsize=xwin, win_ysize=ywin)
+    if maskBand: 
+        data = mb.ReadAsArray(xoff=xoff, yoff=yoff, win_xsize=xwin, win_ysize=ywin)
+    else: 
+        data = sourceBand.ReadAsArray(xoff=xoff, yoff=yoff, win_xsize=xwin, win_ysize=ywin)
+        if dsInfo.scale  != 1.0: data = data*dsInfo.scale
+        if dsInfo.offset != 0.0: data = data+dsInfo.offset
 
     # Correct 'nodata' values
     if autocorrect:
@@ -546,7 +562,7 @@ def isFlipped(source):
     if( dy<0 ): return True
     else: return False
 
-RasterInfo = namedtuple("RasterInfo","srs dtype flipY yAtTop bounds xMin yMin xMax yMax dx dy pixelWidth pixelHeight noData, xWinSize, yWinSize, meta, source")
+RasterInfo = namedtuple("RasterInfo","srs dtype flipY yAtTop bounds xMin yMin xMax yMax dx dy pixelWidth pixelHeight noData, xWinSize, yWinSize, meta, source, scale, offset")
 def rasterInfo(sourceDS):
     """Returns a named tuple containing information relating to the input raster
 
@@ -566,6 +582,8 @@ def rasterInfo(sourceDS):
                     dx:The raster's pixelWidth
                     dy: The raster's pixelHeight
                     noData: The noData value used by the raster
+                    scale: The scale value used by the raster
+                    offset: The offset value used by the raster
                     xWinSize: The width of the raster is pixels
                     yWinSize: The height of the raster is pixels
                     meta: The raster's meta data )
@@ -581,6 +599,8 @@ def rasterInfo(sourceDS):
     sourceBand = sourceDS.GetRasterBand(1)
     output['dtype'] = sourceBand.DataType
     output['noData'] = sourceBand.GetNoDataValue()
+    output['scale'] = sourceBand.GetScale()
+    output['offset'] = sourceBand.GetOffset()
 
     xSize = sourceBand.XSize
     ySize = sourceBand.YSize
@@ -765,11 +785,12 @@ def loadPoint(pt,s):
 
     for xi,yi,ib in zip(xStarts, yStarts, inBounds):
         if not ib:
-            data = np.zeros((window,window))
-            data[:,:] = np.nan
+            data = np.empty((window,window))
         else:
             # Open and read from raster
             data = band.ReadAsArray(xoff=xi, yoff=yi, win_xsize=window, win_ysize=window)
+            if info.scale  != 1.0: data = data*info.scale
+            if info.offset != 0.0: data = data+info.offset
 
             # Look for nodata
             if not info.noData is None:

geokit/core/util.py

@@ -351,7 +351,7 @@ def fitBoundsTo(bounds, dx, dy):
 def quickRaster(bounds, srs, dx, dy, dType="GDT_Byte", noData=None, fill=None, data=None, header=''):
     """GeoKit internal for quickly creating a raster datasource"""
 
-    bounds = fitBoundsTo(bounds, dx, dy)
+    #bounds = fitBoundsTo(bounds, dx, dy)
 
     ## Make a raster dataset and pull the band/maskBand objects
     originX = bounds[0]

geokit/core/vector.py

@@ -962,7 +962,7 @@ def rasterize(source, pixelWidth, pixelHeight, srs=None, bounds=None, where=None
 
         # Do rasterize
         tmp = gdal.Rasterize( output, source, outputBounds=bounds, xRes=pixelWidth, yRes=pixelHeight,
-                              outputSRS=srs, outputType=getattr(gdal, dtype), noData=noData, where=where, 
+                              outputSRS=srs, noData=noData, where=where, 
                               creationOptions=co, targetAlignedPixels=aligned, **kwargs)
 
         if not isRaster(tmp): raise GeoKitRegionMaskError("Rasterization failed!")