diff --git a/bo.txt b/bo.txt new file mode 100644 index 0000000..e69de29 diff --git a/geokit/__init__.py b/geokit/__init__.py index a6261c2..a05c02e 100644 --- a/geokit/__init__.py +++ b/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 diff --git a/geokit/core/extent.py b/geokit/core/extent.py index 6980992..b1719ce 100755 --- a/geokit/core/extent.py +++ b/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) diff --git a/geokit/core/location.py b/geokit/core/location.py index 33e81a7..1d04d9d 100644 --- a/geokit/core/location.py +++ b/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" diff --git a/geokit/core/raster.py b/geokit/core/raster.py index 2531f5e..dbb2614 100644 --- a/geokit/core/raster.py +++ b/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: diff --git a/geokit/core/util.py b/geokit/core/util.py index 05572b1..bd05a86 100644 --- a/geokit/core/util.py +++ b/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] diff --git a/geokit/core/vector.py b/geokit/core/vector.py index 51d1f63..2fa42ab 100644 --- a/geokit/core/vector.py +++ b/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!")