diff --git a/README.md b/README.md index 022254df..6790a106 100755 --- a/README.md +++ b/README.md @@ -144,9 +144,9 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e | [Curvature](xrspatial/curvature.py) | ✅️ | | | ⚠️ | | [Hillshade](xrspatial/hillshade.py) | ✅️ | ✅️ | | | | [Slope](xrspatial/slope.py) | ✅️ | ✅️ | ✅️ | ⚠️ | -| [Terrain Generation](xrspatial/terrain.py) | ✅️ | | | | +| [Terrain Generation](xrspatial/terrain.py) | ✅️ | ✅️ | ✅️ | | | [Viewshed](xrspatial/viewshed.py) | ✅️ | | | | -| [Perlin Noise](xrspatial/perlin.py) | ✅️ | | | | +| [Perlin Noise](xrspatial/perlin.py) | ✅️ | ✅️ | ✅️ | | | [Bump Mapping](xrspatial/bump.py) | ✅️ | | | | ----------- diff --git a/xrspatial/perlin.py b/xrspatial/perlin.py index c6b33325..26a7ea8e 100755 --- a/xrspatial/perlin.py +++ b/xrspatial/perlin.py @@ -1,153 +1,49 @@ -import numpy as np +# std lib +from functools import partial +# 3rd-party +import numpy as np import xarray as xr -from xarray import DataArray - -from xrspatial.utils import ngjit - - -# TODO: change parameters to take agg instead of height / width -def perlin(width: int, - height: int, - freq: tuple = (1, 1), - seed: int = 5) -> xr.DataArray: - """ - Generate perlin noise aggregate. - - Parameters - ---------- - width : int - Width of output aggregate array. - height : int - Height of output aggregate array. - freq : tuple, default=(1,1) - (x, y) frequency multipliers. - seed : int, default=5 - Seed for random number generator. - Returns - ------- - perlin_agg : xarray.DataArray - 2D array of perlin noise values. +try: + import cupy +except ImportError: + class cupy(object): + ndarray = False - References - ---------- - - Paul Panzer: https://stackoverflow.com/questions/42147776/producing-2d-perlin-noise-with-numpy # noqa - - ICA: http://www.mountaincartography.org/mt_hood/pdfs/nighbert_bump1.pdf # noqa - - Examples - -------- - .. plot:: - :include-source: - - import matplotlib.pyplot as plt - from xrspatial import perlin +import dask.array as da +from numba import cuda, jit +import numba as nb - # Generate Perlin Noise Aggregate - perlin_default = perlin(width = 500, height = 300) +# local modules +from xrspatial.utils import has_cuda +from xrspatial.utils import cuda_args - # With Increased x Frequency - perlin_high_x_freq = perlin(width = 500, height = 300, freq = (5, 1)) - # With Increased y Frequency - perlin_high_y_freq = perlin(width = 500, height = 300, freq = (1, 5)) - - # With a Different Seed - perlin_seed_1 = perlin(width = 500, height = 300, seed = 1) - - # Plot Default Perlin - perlin_default.plot(cmap = 'inferno', aspect = 2, size = 4) - plt.title("Default") - - # Plot High x Frequency - perlin_high_x_freq.plot(cmap = 'inferno', aspect = 2, size = 4) - plt.title("High x Frequency") - - # Plot High y Frequency - perlin_high_y_freq.plot(cmap = 'inferno', aspect = 2, size = 4) - plt.title("High y Frequency") - - # Plot Seed = 1 - perlin_seed_1.plot(cmap = 'inferno', aspect = 2, size = 4) - plt.title("Seed = 1") - - .. sourcecode:: python - - >>> print(perlin_default[200:203, 200: 202]) - - array([[0.56800979, 0.56477393], - [0.56651744, 0.56331014], - [0.56499184, 0.56181344]]) - Dimensions without coordinates: y, x - Attributes: - res: 1 - - >>> print(perlin_high_x_freq[200:203, 200: 202]) - - array([[0.5 , 0.48999444], - [0.5 , 0.48999434], - [0.5 , 0.48999425]]) - Dimensions without coordinates: y, x - Attributes: - res: 1 - - >>> print(perlin_high_y_freq[200:203, 200: 202]) - - array([[0.31872961, 0.31756859], - [0.2999256 , 0.2988189 ], - [0.28085118, 0.27979834]]) - Dimensions without coordinates: y, x - Attributes: - res: 1 - - >>> print(perlin_seed_1[200:203, 200: 202]) - - array([[0.12991498, 0.12984185], - [0.13451158, 0.13441514], - [0.13916956, 0.1390495 ]]) - Dimensions without coordinates: y, x - Attributes: - res: 1 - """ - linx = range(width) - liny = range(height) - linx = np.linspace(0, 1, width, endpoint=False) - liny = np.linspace(0, 1, height, endpoint=False) - x, y = np.meshgrid(linx, liny) - data = _perlin(x * freq[0], y * freq[1], seed=seed) - data = (data - np.min(data))/np.ptp(data) - return DataArray(data, dims=['y', 'x'], attrs=dict(res=1)) - - -@ngjit +@jit(nopython=True, nogil=True, parallel=True, cache=True) def _lerp(a, b, x): - return a + x * (b-a) + return a + x * (b - a) -@ngjit +@jit(nopython=True, nogil=True, parallel=True, cache=True) def _fade(t): - return 6 * t**5 - 15 * t**4 + 10 * t**3 + return 6 * t ** 5 - 15 * t ** 4 + 10 * t ** 3 -@ngjit +@jit(nopython=True, nogil=True, parallel=True, cache=True) def _gradient(h, x, y): + # assert(len(h.shape) == 2) vectors = np.array([[0, 1], [0, -1], [1, 0], [-1, 0]]) - dim_ = h.shape - out = np.zeros(dim_) - for j in range(dim_[1]): - for i in range(dim_[0]): + out = np.zeros(h.shape) + for j in nb.prange(h.shape[1]): + for i in nb.prange(h.shape[0]): f = np.mod(h[i, j], 4) g = vectors[f] out[i, j] = g[0] * x[i, j] + g[1] * y[i, j] return out -def _perlin(x, y, seed=0): - np.random.seed(seed) - p = np.arange(2**20, dtype=int) - np.random.shuffle(p) - p = np.stack([p, p]).flatten() - +def _perlin(p, x, y): # coordinates of the top-left xi = x.astype(int) yi = y.astype(int) @@ -161,13 +57,197 @@ def _perlin(x, y, seed=0): v = _fade(yf) # noise components - n00 = _gradient(p[p[xi]+yi], xf, yf) - n01 = _gradient(p[p[xi]+yi+1], xf, yf-1) - n11 = _gradient(p[p[xi+1]+yi+1], xf-1, yf-1) - n10 = _gradient(p[p[xi+1]+yi], xf-1, yf) + n00 = _gradient(p[p[xi] + yi], xf, yf) + n01 = _gradient(p[p[xi] + yi + 1], xf, yf - 1) + n11 = _gradient(p[p[xi + 1] + yi + 1], xf - 1, yf - 1) + n10 = _gradient(p[p[xi + 1] + yi], xf - 1, yf) # combine noises x1 = _lerp(n00, n10, u) x2 = _lerp(n01, n11, u) a = _lerp(x1, x2, v) return a + + +def _perlin_numpy(data: np.ndarray, + freq: tuple, + seed: int) -> np.ndarray: + np.random.seed(seed) + p = np.random.permutation(2**20) + p = np.append(p, p) + + height, width = data.shape + linx = np.linspace(0, freq[0], width, endpoint=False, dtype=np.float32) + liny = np.linspace(0, freq[1], height, endpoint=False, dtype=np.float32) + x, y = np.meshgrid(linx, liny) + + data[:] = _perlin(p, x, y) + data[:] = (data - np.min(data)) / np.ptp(data) + return data + + +def _perlin_dask_numpy(data: da.Array, + freq: tuple, + seed: int) -> da.Array: + np.random.seed(seed) + p = np.random.permutation(2**20) + p = np.append(p, p) + + height, width = data.shape + linx = da.linspace(0, freq[0], width, endpoint=False, dtype=np.float32) + liny = da.linspace(0, freq[1], height, endpoint=False, dtype=np.float32) + x, y = da.meshgrid(linx, liny) + + _func = partial(_perlin, p) + data = da.map_blocks(_func, x, y, meta=np.array((), dtype=np.float32)) + + data = (data - da.min(data)) / da.ptp(data) + return data + + +@cuda.jit(device=True) +def _lerp_gpu(a, b, x): + return a + x * (b - a) + + +@cuda.jit(device=True) +def _fade_gpu(t): + return 6 * t ** 5 - 15 * t ** 4 + 10 * t ** 3 + + +@cuda.jit(device=True) +def _gradient_gpu(vec, h, x, y): + f = h % 4 + return vec[f][0] * x + vec[f][1] * y + + +@cuda.jit(fastmath=True, opt=True) +def _perlin_gpu(p, x0, x1, y0, y1, m, out): + # alloc and initialize array to be used in the gradient routine + vec = cuda.local.array((4, 2), nb.int32) + vec[0][0] = vec[1][0] = vec[2][1] = vec[3][1] = 0 + vec[0][1] = vec[2][0] = 1 + vec[1][1] = vec[3][0] = -1 + + i, j = nb.cuda.grid(2) + if i < out.shape[0] and j < out.shape[1]: + # coordinates of the top-left + y = y0 + i * (y1 - y0) / out.shape[0] + x = x0 + j * (x1 - x0) / out.shape[1] + + # coordinates of the top-left + x_int = int(x) + y_int = int(y) + + # internal coordinates + xf = x - x_int + yf = y - y_int + + # fade factors + u = _fade_gpu(xf) + v = _fade_gpu(yf) + + # noise components + n00 = _gradient_gpu(vec, p[p[x_int] + y_int], xf, yf) + n01 = _gradient_gpu(vec, p[p[x_int] + y_int + 1], xf, yf - 1) + n11 = _gradient_gpu(vec, p[p[x_int + 1] + y_int + 1], xf - 1, yf - 1) + n10 = _gradient_gpu(vec, p[p[x_int + 1] + y_int], xf - 1, yf) + + # combine noises + x1 = _lerp_gpu(n00, n10, u) + x2 = _lerp_gpu(n01, n11, u) + out[i, j] = m * _lerp_gpu(x1, x2, v) + + +def _perlin_cupy(data: cupy.ndarray, + freq: tuple, + seed: int) -> cupy.ndarray: + + # cupy.random.seed(seed) + # p = cupy.random.permutation(2**20) + + # use numpy.random then transfer data to GPU to ensure the same result + # when running numpy backed and cupy backed data array. + np.random.seed(seed) + p = cupy.asarray(np.random.permutation(2**20)) + p = cupy.append(p, p) + + griddim, blockdim = cuda_args(data.shape) + _perlin_gpu[griddim, blockdim](p, 0, freq[0], 0, freq[1], 1, data) + + minimum = cupy.amin(data) + maximum = cupy.amax(data) + data[:] = (data - minimum) / (maximum - minimum) + return data + + +def perlin(agg: xr.DataArray, + freq: tuple = (1, 1), + seed: int = 5, + name: str = 'perlin') -> xr.DataArray: + """ + Generate perlin noise aggregate. + + Parameters + ---------- + agg : xr.DataArray + 2D array of size width x height, will be used to determine + height/ width and which platform to use for calculation. + freq : tuple, default=(1,1) + (x, y) frequency multipliers. + seed : int, default=5 + Seed for random number generator. + + Returns + ------- + perlin_agg : xarray.DataArray + 2D array of perlin noise values. + + References + ---------- + - Paul Panzer: https://stackoverflow.com/questions/42147776/producing-2d-perlin-noise-with-numpy # noqa + - ICA: http://www.mountaincartography.org/mt_hood/pdfs/nighbert_bump1.pdf # noqa + + Examples + -------- + .. plot:: + :include-source: + + import numpy as np + import xarray as xr + from xrspatial import perlin + + W = 4 + H = 3 + data = np.zeros((H, W), dtype=np.float32) + raster = xr.DataArray(data, dims=['y', 'x']) + + perlin_noise = perlin(raster) + + .. sourcecode:: python + + >>> print(perlin_noise) + + array([[0.39268944, 0.27577767, 0.01621884, 0.05518942], + [1. , 0.8229485 , 0.2935367 , 0. ], + [1. , 0.8715414 , 0.41902685, 0.02916668]], dtype=float32) # noqa + Dimensions without coordinates: y, x + """ + + # numpy case + if isinstance(agg.data, np.ndarray): + out = _perlin_numpy(agg.data, freq, seed) + # cupy case + elif has_cuda() and isinstance(agg.data, cupy.ndarray): + out = _perlin_cupy(agg.data, freq, seed) + # dask + numpy case + elif isinstance(agg.data, da.Array): + out = _perlin_dask_numpy(agg.data, freq, seed) + else: + raise TypeError('Unsupported Array Type: {}'.format(type(agg.data))) + + result = xr.DataArray(out, + dims=agg.dims, + attrs=agg.attrs, + name=name) + return result diff --git a/xrspatial/terrain.py b/xrspatial/terrain.py index a38d75ca..657162bd 100755 --- a/xrspatial/terrain.py +++ b/xrspatial/terrain.py @@ -1,22 +1,189 @@ +# std lib +from functools import partial +from typing import Union, Tuple, List +from typing import Optional + +# 3rd-party import numpy as np +import xarray as xr import pandas as pd import datashader as ds -from typing import Optional -import xarray as xr -from xarray import DataArray -from .perlin import _perlin +try: + import cupy +except ImportError: + class cupy(object): + ndarray = False + +import dask.array as da + +# local modules +from xrspatial.utils import has_cuda +from xrspatial.utils import cuda_args +from xrspatial.utils import get_dataarray_resolution +from .perlin import _perlin, _perlin_gpu + + +def _scale(value, old_range, new_range): + d = (value - old_range[0]) / (old_range[1] - old_range[0]) + return d * (new_range[1] - new_range[0]) + new_range[0] + + +def _gen_terrain(height_map, seed, x_range=(0, 1), y_range=(0, 1)): + height, width = height_map.shape + + # multiplier, (xfreq, yfreq) + NOISE_LAYERS = ((1 / 2**i, (2**i, 2**i)) for i in range(16)) + + linx = np.linspace( + x_range[0], x_range[1], width, endpoint=False, dtype=np.float32 + ) + liny = np.linspace( + y_range[0], y_range[1], height, endpoint=False, dtype=np.float32 + ) + x, y = np.meshgrid(linx, liny) + nrange = np.arange(2**20, dtype=int) + for i, (m, (xfreq, yfreq)) in enumerate(NOISE_LAYERS): + np.random.seed(seed+i) + p = np.random.permutation(nrange) + p = np.append(p, p) + + noise = _perlin(p, x * xfreq, y * yfreq) * m + + height_map += noise + + height_map /= (1.00 + 0.50 + 0.25 + 0.13 + 0.06 + 0.03) + height_map = height_map ** 3 + return height_map + + +def _terrain_numpy(data: np.ndarray, + seed: int, + x_range_scaled: tuple, + y_range_scaled: tuple, + zfactor: int) -> np.ndarray: + + # data.fill(0) + data = data * 0 + data[:] = _gen_terrain( + data, seed, x_range=x_range_scaled, y_range=y_range_scaled + ) + + data = (data - np.min(data))/np.ptp(data) + data[data < 0.3] = 0 # create water + data *= zfactor + + return data + + +def _terrain_dask_numpy(data: da.Array, + seed: int, + x_range_scaled: tuple, + y_range_scaled: tuple, + zfactor: int) -> da.Array: + data = data * 0 + + height, width = data.shape + linx = da.linspace( + x_range_scaled[0], x_range_scaled[1], width, endpoint=False, + dtype=np.float32 + ) + liny = da.linspace( + y_range_scaled[0], y_range_scaled[1], height, endpoint=False, + dtype=np.float32 + ) + x, y = da.meshgrid(linx, liny) + + nrange = np.arange(2 ** 20, dtype=int) + + # multiplier, (xfreq, yfreq) + NOISE_LAYERS = ((1 / 2 ** i, (2 ** i, 2 ** i)) for i in range(16)) + for i, (m, (xfreq, yfreq)) in enumerate(NOISE_LAYERS): + np.random.seed(seed + i) + p = np.random.permutation(nrange) + p = np.append(p, p) + + _func = partial(_perlin, p) + noise = da.map_blocks( + _func, + x * xfreq, + y * yfreq, + meta=np.array((), dtype=np.float32) + ) + + data += noise * m + + data /= (1.00 + 0.50 + 0.25 + 0.13 + 0.06 + 0.03) + data = data ** 3 + + data = (data - np.min(data)) / np.ptp(data) + data[data < 0.3] = 0 # create water + data *= zfactor + + return data + + +def _terrain_gpu(height_map, seed, x_range=(0, 1), y_range=(0, 1)): + + NOISE_LAYERS = ((1 / 2**i, (2**i, 2**i)) for i in range(16)) + + noise = cupy.empty_like(height_map, dtype=np.float32) + + griddim, blockdim = cuda_args(height_map.shape) + + for i, (m, (xfreq, yfreq)) in enumerate(NOISE_LAYERS): + + # cupy.random.seed(seed+i) + # p = cupy.random.permutation(2**20) + + # use numpy.random then transfer data to GPU to ensure the same result + # when running numpy backed and cupy backed data array. + np.random.seed(seed+i) + p = cupy.asarray(np.random.permutation(2**20)) + p = cupy.append(p, p) + + _perlin_gpu[griddim, blockdim]( + p, x_range[0] * xfreq, x_range[1] * xfreq, + y_range[0] * yfreq, y_range[1] * yfreq, + m, noise + ) + + height_map += noise + + height_map /= (1.00 + 0.50 + 0.25 + 0.13 + 0.06 + 0.03) + height_map = height_map ** 3 + return height_map + + +def _terrain_cupy(data: cupy.ndarray, + seed: int, + x_range_scaled: tuple, + y_range_scaled: tuple, + zfactor: int) -> cupy.ndarray: + + data = data * 0 + + data[:] = _terrain_gpu(data, seed, + x_range=x_range_scaled, + y_range=y_range_scaled) + minimum = cupy.amin(data) + maximum = cupy.amax(data) + data[:] = (data - minimum)/(maximum - minimum) + data[data < 0.3] = 0 # create water + data *= zfactor + + return data -# TODO: add optional name parameter `name='terrain'` -def generate_terrain(x_range: tuple = (0, 500), + +def generate_terrain(agg: xr.DataArray, + x_range: tuple = (0, 500), y_range: tuple = (0, 500), - width: int = 25, - height: int = 30, - canvas: ds.Canvas = None, seed: int = 10, zfactor: int = 4000, - full_extent: Optional[str] = None) -> xr.DataArray: + full_extent: Optional[Union[Tuple, List]] = None, + name: str = 'terrain' + ) -> xr.DataArray: """ Generates a pseudo-random terrain which can be helpful for testing raster functions. @@ -27,12 +194,6 @@ def generate_terrain(x_range: tuple = (0, 500), Range of x values. x_range : tuple, default=(0, 500) Range of y values. - width : int, default=25 - Width of output data array in pixels. - height : int, default=30 - Height of output data array in pixels. - canvas : ds.Canvas, default=None - Instance for passing output dimensions / ranges. seed : int, default=10 Seed for random number generator. zfactor : int, default=4000 @@ -55,133 +216,83 @@ def generate_terrain(x_range: tuple = (0, 500), .. plot:: :include-source: - import datashader as ds - import matplotlib.pyplot as plt - from xrspatial import generate_terrain, aspect - - # Create Canvas - W = 500 - H = 300 - cvs = ds.Canvas(plot_width = W, - plot_height = H, - x_range = (-20e6, 20e6), - y_range = (-20e6, 20e6)) - - # Generate Example Terrain - terrain_agg = generate_terrain(canvas = cvs) - - # Edit Attributes - terrain_agg = terrain_agg.assign_attrs( - { - 'Description': 'Example Terrain', - 'units': 'km', - 'Max Elevation': '4000', - } - ) + import numpy as np + import xarray as xr + from xrspatial import generate_terrain + + W = 4 + H = 3 + data = np.zeros((H, W), dtype=np.float32) + raster = xr.DataArray(data, dims=['y', 'x']) - terrain_agg = terrain_agg.rename({'x': 'lon', 'y': 'lat'}) - terrain_agg = terrain_agg.rename('Elevation') + xrange = (-20e6, 20e6) + yrange = (-20e6, 20e6) + seed = 2 + zfactor = 10 - # Plot Terrain - terrain_agg.plot(cmap = 'terrain', aspect = 2, size = 4) - plt.title("Terrain") - plt.ylabel("latitude") - plt.xlabel("longitude") + terrain = generate_terrain(raster, xrange, yrange, seed, zfactor) .. sourcecode:: python - >>> print(terrain_agg[200:203, 200:202]) - - array([[1264.02249454, 1261.94748873], - [1285.37061171, 1282.48046696], - [1306.02305679, 1303.40657515]]) + >>> print(terrain) + + array([[ 6.8067746, 5.263137 , 4.664292 , 6.821344 ], + [ 7.4834156, 6.9849734, 4.3545456, 0. ], + [ 6.7674546, 10. , 7.0946655, 7.015267 ]], dtype=float32) # noqa Coordinates: - * lon (lon) float64 -3.96e+06 -3.88e+06 - * lat (lat) float64 6.733e+06 6.867e+06 7e+06 + * x (x) float64 -1.5e+07 -5e+06 5e+06 1.5e+07 + * y (y) float64 -1.333e+07 0.0 1.333e+07 Attributes: - res: 1 - Description: Example Terrain - units: km - Max Elevation: 4000 + res: (10000000.0, -13333333.333333334) """ - def _gen_heights(bumps): - out = np.zeros(len(bumps)) - for i, b in enumerate(bumps): - x = b[0] - y = b[1] - val = agg.data[y, x] - if val >= 0.33 and val <= 3: - out[i] = 0.1 - return out - - def _scale(value, old_range, new_range): - d = (value - old_range[0]) / (old_range[1] - old_range[0]) - return d * (new_range[1] - new_range[0]) + new_range[0] - - mercator_extent = (-np.pi * 6378137, -np.pi * 6378137, - np.pi * 6378137, np.pi * 6378137) - crs_extents = {'3857': mercator_extent} - - if isinstance(full_extent, str): - full_extent = crs_extents[full_extent] - - elif full_extent is None: - full_extent = (canvas.x_range[0], canvas.y_range[0], - canvas.x_range[1], canvas.y_range[1]) + + height, width = agg.shape + + if full_extent is None: + full_extent = (x_range[0], y_range[0], + x_range[1], y_range[1]) elif not isinstance(full_extent, (list, tuple)) and len(full_extent) != 4: - raise TypeError('full_extent must be tuple(4) or str wkid') + raise TypeError('full_extent must be tuple(4)') full_xrange = (full_extent[0], full_extent[2]) full_yrange = (full_extent[1], full_extent[3]) - x_range_scaled = (_scale(canvas.x_range[0], full_xrange, (0.0, 1.0)), - _scale(canvas.x_range[1], full_xrange, (0.0, 1.0))) + x_range_scaled = (_scale(x_range[0], full_xrange, (0.0, 1.0)), + _scale(x_range[1], full_xrange, (0.0, 1.0))) - y_range_scaled = (_scale(canvas.y_range[0], full_yrange, (0.0, 1.0)), - _scale(canvas.y_range[1], full_yrange, (0.0, 1.0))) + y_range_scaled = (_scale(y_range[0], full_yrange, (0.0, 1.0)), + _scale(y_range[1], full_yrange, (0.0, 1.0))) - data = _gen_terrain(canvas.plot_width, canvas.plot_height, seed, - x_range=x_range_scaled, y_range=y_range_scaled) + # numpy case + if isinstance(agg.data, np.ndarray): + out = _terrain_numpy( + agg.data, seed, x_range_scaled, y_range_scaled, zfactor + ) + # cupy case + elif has_cuda() and isinstance(agg.data, cupy.ndarray): + out = _terrain_cupy( + agg.data, seed, x_range_scaled, y_range_scaled, zfactor + ) + # dask + numpy case + elif isinstance(agg.data, da.Array): + out = _terrain_dask_numpy( + agg.data, seed, x_range_scaled, y_range_scaled, zfactor + ) + else: + raise TypeError('Unsupported Array Type: {}'.format(type(agg.data))) - data = (data - np.min(data))/np.ptp(data) - data[data < 0.3] = 0 # create water - data *= zfactor + canvas = ds.Canvas( + plot_width=width, plot_height=height, x_range=x_range, y_range=y_range + ) # DataArray coords were coming back different from cvs.points... hack_agg = canvas.points(pd.DataFrame({'x': [], 'y': []}), 'x', 'y') - agg = DataArray(data, - name='terrain', - coords=hack_agg.coords, - dims=hack_agg.dims, - attrs={'res': 1}) - - return agg - - -def _gen_terrain(width, height, seed, x_range=None, y_range=None): - - if not x_range: - x_range = (0, 1) - - if not y_range: - y_range = (0, 1) - - # multiplier, (xfreq, yfreq) - NOISE_LAYERS = ((1 / 2**i, (2**i, 2**i)) for i in range(16)) - - linx = np.linspace(x_range[0], x_range[1], width, endpoint=False) - liny = np.linspace(y_range[0], y_range[1], height, endpoint=False) - x, y = np.meshgrid(linx, liny) - - height_map = None - for i, (m, (xfreq, yfreq)) in enumerate(NOISE_LAYERS): - noise = _perlin(x * xfreq, y * yfreq, seed=seed + i) * m - if height_map is None: - height_map = noise - else: - height_map += noise - - height_map /= (1.00 + 0.50 + 0.25 + 0.13 + 0.06 + 0.03) - height_map = height_map ** 3 - return height_map + res = get_dataarray_resolution(hack_agg) + result = xr.DataArray(out, + name=name, + coords=hack_agg.coords, + dims=hack_agg.dims, + attrs={'res': res}) + + return result diff --git a/xrspatial/tests/test_perlin.py b/xrspatial/tests/test_perlin.py old mode 100644 new mode 100755 index a15623ad..eabb3fca --- a/xrspatial/tests/test_perlin.py +++ b/xrspatial/tests/test_perlin.py @@ -1,6 +1,58 @@ +import pytest +import numpy as np +import xarray as xr +import dask.array as da + +from xrspatial.utils import has_cuda +from xrspatial.utils import doesnt_have_cuda from xrspatial import perlin -def test_perlin(): - print('Perlin ({}) test not implemented...'.format(perlin)) +def create_test_arr(backend='numpy'): + W = 50 + H = 50 + data = np.zeros((H, W), dtype=np.float32) + raster = xr.DataArray(data, dims=['y', 'x']) + + if has_cuda() and 'cupy' in backend: + import cupy + raster.data = cupy.asarray(raster.data) + + if 'dask' in backend: + raster.data = da.from_array(raster.data, chunks=(10, 10)) + + return raster + + +def test_perlin_cpu(): + # vanilla numpy version + data_numpy = create_test_arr() + perlin_numpy = perlin(data_numpy) + + # dask + data_dask = create_test_arr(backend='dask') + perlin_dask = perlin(data_dask) + assert isinstance(perlin_dask.data, da.Array) + + perlin_dask = perlin_dask.compute() + assert np.isclose( + perlin_numpy.data, perlin_dask.data, + rtol=1e-05, atol=1e-07, equal_nan=True + ).all() + + +@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available") +def test_perlin_gpu(): + # vanilla numpy version + data_numpy = create_test_arr() + perlin_numpy = perlin(data_numpy) + + # cupy + data_cupy = create_test_arr(backend='cupy') + perlin_cupy = perlin(data_cupy) + + assert np.isclose( + perlin_numpy.data, perlin_cupy.data, + rtol=1e-05, atol=1e-07, equal_nan=True + ).all() diff --git a/xrspatial/tests/test_terrain.py b/xrspatial/tests/test_terrain.py old mode 100644 new mode 100755 index f1389a39..476d5611 --- a/xrspatial/tests/test_terrain.py +++ b/xrspatial/tests/test_terrain.py @@ -1,16 +1,58 @@ -import datashader as ds +import pytest +import numpy as np +import xarray as xr +import dask.array as da + +from xrspatial.utils import has_cuda +from xrspatial.utils import doesnt_have_cuda from xrspatial import generate_terrain -W = 25 -H = 30 -X_RANGE = (0, 500) -Y_RANGE = (0, 500) +def create_test_arr(backend='numpy'): + W = 50 + H = 50 + data = np.zeros((H, W), dtype=np.float32) + raster = xr.DataArray(data, dims=['y', 'x']) + + if has_cuda() and 'cupy' in backend: + import cupy + raster.data = cupy.asarray(raster.data) + + if 'dask' in backend: + raster.data = da.from_array(raster.data, chunks=(10, 10)) + + return raster + + +def test_terrain_cpu(): + # vanilla numpy version + data_numpy = create_test_arr() + terrain_numpy = generate_terrain(data_numpy) + + # dask + data_dask = create_test_arr(backend='dask') + terrain_dask = generate_terrain(data_dask) + assert isinstance(terrain_dask.data, da.Array) + + terrain_dask = terrain_dask.compute() + assert np.isclose( + terrain_numpy.data, terrain_dask.data, + rtol=1e-05, atol=1e-07, equal_nan=True + ).all() + + +@pytest.mark.skipif(doesnt_have_cuda(), reason="CUDA Device not Available") +def test_terrain_gpu(): + # vanilla numpy version + data_numpy = create_test_arr() + terrain_numpy = generate_terrain(data_numpy) + # cupy + data_cupy = create_test_arr(backend='cupy') + terrain_cupy = generate_terrain(data_cupy) -def test_generate_terrain(): - csv = ds.Canvas(x_range=X_RANGE, y_range=Y_RANGE, - plot_width=W, plot_height=H) - terrain = generate_terrain(canvas=csv) - assert terrain is not None + assert np.isclose( + terrain_numpy.data, terrain_cupy.data, + rtol=1e-05, atol=1e-07, equal_nan=True + ).all()