From 7da588573c8c0ada6ab8f874f40bd77841f45936 Mon Sep 17 00:00:00 2001 From: Chris B Date: Mon, 10 Sep 2018 13:45:02 +0100 Subject: [PATCH] Revert "WIP: Geo-related transfer functions (#624)" (#652) This reverts commit 40f61f8406527ce23e352eae338f71c183e03914. (#624 was merged prematurely.) --- datashader/geo.py | 311 -------------- datashader/tests/test_geo.py | 157 ------- examples/user_guide/8_Geography.ipynb | 596 +------------------------- tox.ini | 1 - 4 files changed, 2 insertions(+), 1063 deletions(-) delete mode 100644 datashader/geo.py delete mode 100644 datashader/tests/test_geo.py diff --git a/datashader/geo.py b/datashader/geo.py deleted file mode 100644 index 29874e2af..000000000 --- a/datashader/geo.py +++ /dev/null @@ -1,311 +0,0 @@ -""" -This module contains geoscience-related transfer functions whose use is completely optional. -""" - -from __future__ import division - -import numpy as np -import datashader.transfer_functions as tf - -from datashader.colors import rgb -from datashader.utils import ngjit -from xarray import DataArray - -__all__ = ['mean', 'binary', 'slope', 'aspect', 'ndvi', 'hillshade'] - - -def _shade(altituderad, aspect, azimuthrad, slope): - shade = np.sin(altituderad) * np.sin(slope) + np.cos(altituderad) * np.cos(slope) * np.cos(azimuthrad - aspect) - return shade - -def _threshold_hs(img): - dt = np.dtype((np.int32, {'r': (np.uint8, 0), - 'g': (np.uint8, 1), - 'b': (np.uint8, 2), - 'a': (np.uint8, 3)})) - img.data = np.where(img.data.view(dtype=dt)['r'] > 105, 0, img) # awk. - return img - -def _simple_hs(altitude, aspect, azimuth, slope, cmap, alpha, out_type='image'): - _shaded = _shade(altitude, aspect, azimuth, slope) - - agg = DataArray(_shaded, dims=['y','x']) - agg.data = np.where(agg > agg.mean(), 0, agg) - if out_type == 'image': - img = tf.shade(agg, cmap=cmap, how='linear', alpha=alpha) - return _threshold_hs(img) - elif out_type == 'data': - return agg - else: - raise ValueError("Unknown out_type: {0}".format(out_type)) - - -def _mdow_hs(altitude, aspect, azimuth, slope, cmap, alpha, out_type='image'): - alt = np.deg2rad(30) - shade = np.sum([_shade(alt, aspect, np.deg2rad(225), slope), - _shade(alt, aspect, np.deg2rad(270), slope), - _shade(alt, aspect, np.deg2rad(315), slope), - _shade(alt, aspect, np.deg2rad(360), slope)], axis=0) - shade /= 4 - - agg = DataArray(shade, dims=['y', 'x']) - agg.data = np.where(agg > agg.mean(), 0, agg) - if out_type == 'image': - img = tf.shade(agg, cmap=cmap, how='linear', alpha=alpha) - return _threshold_hs(img) - elif out_type == 'data': - return agg - else: - raise ValueError("Unknown out_type: {0}".format(out_type)) - - -_hillshade_lookup = {'simple': _simple_hs, - 'mdow': _mdow_hs} - -def _normalize_hillshade_how(how): - if how in _hillshade_lookup: - return _hillshade_lookup[how] - raise ValueError("Unknown hillshade method: {0}".format(how)) - -def hillshade(agg, - altitude=30, - azimuth=315, - alpha=70, - how='mdow', - out_type='image', - cmap=['#C7C7C7', '#000000']): - """Convert a 2D DataArray to an hillshaded image with specified colormap. - - Parameters - ---------- - agg : DataArray - altitude : int, optional (default: 30) - Altitude angle of the sun specified in degrees. - azimuth : int, optional (default: 315) - The angle between the north vector and the perpendicular projection - of the light source down onto the horizon specified in degrees. - cmap : list of colors or matplotlib.colors.Colormap, optional - The colormap to use. Can be either a list of colors (in any of the - formats described above), or a matplotlib colormap object. - Default is `["lightgray", "black"]` - how : str or callable, optional - The hillshade method to use. Valid strings are 'mdow' [default], - 'simple'. - alpha : int, optional - Value between 0 - 255 representing the alpha value of pixels which contain - data (i.e. non-nan values). Regardless of this value, `NaN` values are - set to fully transparent. - - Returns - ------- - Datashader Image - - """ - global _threshold_hs - - if not isinstance(agg, DataArray): - raise TypeError("agg must be instance of DataArray") - - azimuthrad = np.deg2rad(azimuth) - altituderad = np.deg2rad(altitude) - y, x = np.gradient(agg.data) - slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y)) - aspect = np.arctan2(-y, x) - how = _normalize_hillshade_how(how) - return how(altituderad, aspect, azimuthrad, slope, cmap, alpha, out_type) - - -@ngjit -def _horn_slope(data, cellsize, use_percent=True): - out = np.zeros_like(data) - rows, cols = data.shape - if use_percent: - for y in range(1, rows-1): - for x in range(1, cols-1): - a,b,c,d,f,g,h,i = [data[y-1, x-1], data[y, x-1], data[y+1, x-1], - data[y-1, x], data[y+1, x], - data[y-1, x+1], data[y, x+1], data[y+1, x+1]] - dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize) - dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize) - out[y, x] = (dz_dx ** 2 + dz_dy ** 2) ** .5 - else: - for y in range(1, rows-1): - for x in range(1, cols-1): - a,b,c,d,f,g,h,i = [data[y-1, x-1], data[y, x-1], data[y+1, x-1], - data[y-1, x], data[y+1, x], - data[y-1, x+1], data[y, x+1], data[y+1, x+1]] #NOQA - dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize) - dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize) - p = (dz_dx ** 2 + dz_dy ** 2) ** .5 - out[y, x] = np.arctan(p) * 180 / np.pi - - return out - - -def slope(agg, units='percent'): - """Returns slope of input aggregate in percent. - - Parameters - ---------- - agg : DataArray - units : str, optional (default: percent) - The units of the return values. options `percent`, or `degrees`. - - Returns - ------- - data: DataArray - """ - - if not isinstance(agg, DataArray): - raise TypeError("agg must be instance of DataArray") - - if units not in ('percent', 'degree'): - raise ValueError('Invalid slope units: options (percent, degree)') - - if not agg.attrs.get('res'): - #TODO: maybe monkey-patch a "res" attribute valueing unity is reasonable - raise ValueError('input xarray must have numeric `res` attr.') - - use_percent = units == 'percent' - slope_agg = _horn_slope(agg.data, - agg.attrs['res'], - use_percent=use_percent) - - return DataArray(slope_agg, - name='slope', - coords=agg.coords, - dims=agg.dims, - attrs=agg.attrs) - -@ngjit -def _ndvi(nir_data, red_data): - out = np.zeros_like(nir_data) - rows, cols = nir_data.shape - for y in range(0, rows): - for x in range(0, cols): - nir = nir_data[y, x] - red = red_data[y, x] - soma = nir + red - if soma != 0: - out[y, x] = (nir - red) / soma - return out - -def ndvi(nir_agg, red_agg, units='percent'): - """Returns Normalized Difference Vegetation Index (NDVI). - - Parameters - ---------- - nir_agg : DataArray - near-infrared band data - red_agg : DataArray - red band data - - Returns - ------- - data: DataArray - """ - - if not isinstance(nir_agg, DataArray): - raise TypeError("nir_agg must be instance of DataArray") - - if not isinstance(red_agg, DataArray): - raise TypeError("red_agg must be instance of DataArray") - - if not red_agg.shape == nir_agg.shape: - raise ValueError("red_agg and nir_agg expected to have equal shapes") - - return DataArray(_ndvi(nir_agg.data, red_agg.data), - attrs=nir_agg.attrs) - -@ngjit -def _horn_aspect(data): - out = np.zeros_like(data) - rows, cols = data.shape - for y in range(1, rows-1): - for x in range(1, cols-1): - a,b,c,d,f,g,h,i = [data[y-1, x-1], data[y, x-1], data[y+1, x-1], - data[y-1, x], data[y+1, x], - data[y-1, x+1], data[y, x+1], data[y+1, x+1]] - dz_dx = ((c + 2 * f + i) - (a + 2 * d + g)) - dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) - aspect = np.arctan2(dz_dy, -dz_dx) * 180 / np.pi - out[y, x] = aspect + 180 - return out - - -def aspect(agg): - """Returns downward slope direction in compass degrees (0 - 360) with 0 at 12 o'clock. - - Parameters - ---------- - agg : DataArray - - Returns - ------- - data: DataArray - """ - - if not isinstance(agg, DataArray): - raise TypeError("agg must be instance of DataArray") - - return DataArray(_horn_aspect(agg.data), - dims=['y', 'x'], - attrs=agg.attrs) - - -def color_values(agg, color_key, alpha=255): - - def _convert_color(c): - r, g, b = rgb(c) - return np.array([r, g, b, alpha]).astype(np.uint8).view(np.uint32)[0] - - _converted_colors = {k: _convert_color(v) for k, v in color_key.items()} - f = np.vectorize(lambda v: _converted_colors.get(v, 0)) - return tf.Image(f(agg.data)) - - -@ngjit -def _binary(data, values): - out = np.zeros_like(data) - rows, cols = data.shape - for x in range(0, rows): - for y in range(0, cols): - if data[y, x] in values: - out[y, x] = True - else: - out[y, x] = False - return out - -def binary(agg, values): - return DataArray(_binary(agg.data, values), - dims=['y', 'x'], - attrs=agg.attrs) - - -@ngjit -def _mean(data): - out = np.zeros_like(data) - rows, cols = data.shape - for y in range(1, rows-1): - for x in range(1, cols-1): - a,b,c,d,e,f,g,h,i = [data[y-1, x-1], data[y, x-1], data[y+1, x-1], - data[y-1, x], data[y, x], data[y+1, x], - data[y-1, x+1], data[y, x+1], data[y+1, x+1]] - - out[y, x] = (a+b+c+d+e+f+g+h+i) / 9 - return out - - -def mean(agg, pad=None): - """ - Returns Mean filtered array using a 3x3 window - - Parameters - ---------- - agg : DataArray - - Returns - ------- - data: DataArray - """ - return DataArray(_mean(agg.data), dims=['y', 'x'], attrs=agg.attrs) diff --git a/datashader/tests/test_geo.py b/datashader/tests/test_geo.py deleted file mode 100644 index c22d7d025..000000000 --- a/datashader/tests/test_geo.py +++ /dev/null @@ -1,157 +0,0 @@ -# import pytest - -from os import path - -import datashader as ds -import xarray as xr -import numpy as np - -from datashader import geo #mean, binary, slope, aspect, ndvi, hillshade - -# use landsat data from sample data -BASE_PATH = path.split(__file__)[0] -DATA_PATH = path.abspath(path.join(BASE_PATH, 'data')) -# TEST_RASTER_PATH = path.join(DATA_PATH, 'landsat.tif???') # TODO: update me! -TEST_RASTER_PATH = path.join(DATA_PATH, 'world.rgb.tif') - - -# ----- -# Utils -# -with xr.open_rasterio(TEST_RASTER_PATH) as src: - res = ds.utils.calc_res(src) - left, bottom, right, top = ds.utils.calc_bbox(src.x.values, src.y.values, res) - cvs = ds.Canvas(plot_width=2, - plot_height=2, - x_range=(left, right), - y_range=(bottom, top)) - -def _raster_aggregate_default(): - with xr.open_rasterio(TEST_RASTER_PATH) as src: - agg = cvs.raster(src) - assert agg is not None - return agg - -# Take a numpy array and make it randomly sparse (50% -> zeros) -# -def _do_sparse_array(data_array): - import random - indx = list(zip(*np.where(data_array))) - pos = random.sample(range(data_array.size), data_array.size//2) - indx = np.asarray(indx)[pos] - r = indx[:,0] - c = indx[:,1] - data_half = data_array.copy() - data_half[r,c] = 0 - return data_half - -def _do_gaussian_array(): - _x = np.linspace(0, 50, 101) - _y = _x.copy() - _mean = 25 - _sdev = 5 - X,Y = np.meshgrid(_x, _y, sparse=True) - x_fac = -np.power(X-_mean, 2) - y_fac = -np.power(Y-_mean, 2) - gaussian = np.exp((x_fac+y_fac)/(2*_sdev**2)) / (2.5*_sdev) - return gaussian -# -# ----- - -data_random = np.random.random_sample((100,100)) -data_random_sparse = _do_sparse_array(data_random) -data_gaussian = _do_gaussian_array() - - -def test_mean_transfer_function(): - da = xr.DataArray(data_random) - da_mean = geo.mean(da) - assert da.shape == da_mean.shape - - # Overall mean value should be the same as the original array. - # Considering the default behaviour to 'mean' is to pad the borders - # with zeros, the mean value of the filtered array will be slightly - # smaller (considering 'data_random' is positive). - assert da_mean.mean() <= data_random.mean() - - # And if we pad the borders with the original values, we should have a - # 'mean' filtered array with _mean_ value very similar to the original one. - da_mean[0,:] = data_random[0,:] - da_mean[-1,:]= data_random[-1,:] - da_mean[:,0] = data_random[:,0] - da_mean[:,-1]= data_random[:,-1] - assert abs(da_mean.mean() - data_random.mean()) < 10**-3 - -def test_slope_transfer_function(): - """ - Assert slope transfer function - """ - da = xr.DataArray(data_gaussian, attrs={'res':1}) - da_slope = geo.slope(da) - assert da.shape == da_slope.shape - - assert da_slope.sum() > 0 - - # In the middle of the array, there is the maximum of the gaussian; - # And there the slope must be zero. - _imax = np.where(da == da.max()) - assert da_slope[_imax] == 0 - -def test_aspect_transfer_function(): - """ - Assert aspect transfer function - """ - da = xr.DataArray(data_gaussian, attrs={'res':1}) - da_aspect = geo.aspect(da) - assert da.shape == da_aspect.shape - - # Running clockwise, from [0:360] degrees, with origin the the vertical axis - y_mid = data_gaussian.shape[0]//2 - x_mid = data_gaussian.shape[1]//2 - - print(da_aspect) - # middle-top - assert da_aspect[1,x_mid] == 0 or da_aspect[1,x_mid] == 360 - # right-middle - assert da_aspect[y_mid,-2] == 90 - # middle-bottom - assert da_aspect[-2,x_mid] == 180 - # left-middle - assert da_aspect[y_mid,1] == 270 - # top-right - assert da_aspect[1,-2] == 45 - # bottom-right - assert da_aspect[-2,-2] == 135 - # bottom-left - assert da_aspect[-2,1] == 225 - # top-left - assert da_aspect[1,1] == 315 - -def test_hillshade_simple_transfer_function(): - """ - Assert Simple Hillshade transfer function - """ - da_gaussian = xr.DataArray(data_gaussian) - da_gaussian_shade = geo.hillshade(da_gaussian, how='simple', out_type='data') - - assert da_gaussian_shade.mean() > 0 - assert da_gaussian_shade[1,1] == 0 - assert da_gaussian_shade[60,60] > 0 - -def test_ndvi_transfer_function(): - """ - Assert aspect transfer function - """ - _x = np.mgrid[1:0:21j] - a,b = np.meshgrid(_x,_x) - red = a*b - nir = (a*b)[::-1,::-1] - - da_red = xr.DataArray(red) - da_nir = xr.DataArray(nir) - da_ndvi = geo.ndvi(da_nir, da_red) - - assert da_ndvi[0,0] == -1 - assert da_ndvi[-1,-1] == 1 - assert da_ndvi[5,10] == da_ndvi[10,5] == -0.5 - assert da_ndvi[15,10] == da_ndvi[10,15] == 0.5 diff --git a/examples/user_guide/8_Geography.ipynb b/examples/user_guide/8_Geography.ipynb index f1ab4ad80..c1ee5a594 100644 --- a/examples/user_guide/8_Geography.ipynb +++ b/examples/user_guide/8_Geography.ipynb @@ -4,606 +4,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# 8. Geography\n", - "\n", - "The functions are part of `datashader.geo` module.\n", - "\n", - "* [Slope](#ds.geo---slope-function)\n", - "* [Aspect](#ds.geo---aspect-function)\n", - "* [Hillshade](#ds.geo---hillshade-function)\n", - "* [NDVI](#ds.geo---ndvi-function)\n", - "* [Mean](#ds.geo---mean-function)" + "This notebook is under construction; please see [geoviews.org](http://geoviews.org) for several examples of using datashader with geoographical data." ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from datashader.geo import slope, aspect, hillshade, ndvi, binary" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import datashader as ds\n", - "import xarray as xr\n", - "\n", - "import holoviews as hv\n", - "hv.extension('matplotlib','bokeh')\n", - "COLORMAP = 'copper'\n", - "\n", - "from IPython.display import display, Markdown\n", - "\n", - "def version(mod):\n", - " print('{} version: {}'.format(mod.__name__,mod.__version__))\n", - "\n", - "version(ds)\n", - "version(xr)\n", - "version(hv)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def helper_normalize(data_array):\n", - " \"\"\"\n", - " Simple normalization helper\n", - " \"\"\"\n", - " return (data_array - data_array.min()) / (data_array.max() - data_array.min())\n", - "\n", - "\n", - "def helper_do_random():\n", - " \"\"\"Create a uniformly random [0:1] array\"\"\"\n", - " import numpy as np\n", - " import xarray as xr\n", - "\n", - " data_random = np.random.random_sample((Y_SIZE,X_SIZE))\n", - " \n", - " return xr.DataArray(data_random, attrs={'res':1})\n", - "\n", - "\n", - "def helper_do_gaussian():\n", - " \"\"\"Create a gaussian surface\"\"\"\n", - " import numpy as np\n", - " import xarray as xr\n", - "\n", - " X = Y = np.linspace(-4, 4, 301)\n", - " X,Y = np.meshgrid(X, Y, sparse=True)\n", - " X_fac = -np.power(X, 2)\n", - " Y_fac = -np.power(Y, 2)\n", - "\n", - " data_gaussian = np.exp((X_fac+Y_fac))\n", - " data_gaussian = helper_normalize(data_gaussian)\n", - "\n", - " return xr.DataArray(data_gaussian, attrs={'res':1})\n", - "\n", - "\n", - "def helper_do_gaussian_inverse():\n", - " \"\"\"Return a gaussian with inverted values\"\"\"\n", - " data_gaussian = helper_do_gaussian()\n", - " data_gaussian_inv = -1 * (data_gaussian - 1)\n", - " return data_gaussian_inv" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Slope\n", - "* Burrough, P. A., and McDonell, R. A., 1998. *Principles of Geographical Information Systems* (Oxford University Press, New York), pp 406\n", - "\n", - "[Slope](https://en.wikipedia.org/wiki/Slope) is the inclination of a surface. \n", - "In geography, *slope* is amount of change in elevation of a terrain regarding its surroundings." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "display(Markdown(slope.__doc__))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# We will use a gaussian surface to show the output of the `slope` function.\n", - "\n", - "xarray_gaussian = helper_do_gaussian()\n", - "\n", - "da_slope = slope(xarray_gaussian)\n", - "\n", - "img_gauss = hv.Image(xarray_gaussian.data, label='Gaussian curve').options(colorbar=True).options(cmap=COLORMAP)\n", - "img_slope = hv.Image(da_slope.data, label='Slope output').options(colorbar=True).options(cmap=COLORMAP)\n", - "img_slope + img_gauss" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Data\n", - "\n", - "The following section provides synthetic and actual geospatial data to demonstrate geoscience functionality.\n", - "\n", - "\n", - "#### Synthetic\n", - "\n", - "* [Gaussian](#Gaussian) : a gaussian surface that will be used for the `geo.hillshade` function\n", - " * Since a Gaussian distribution is simple, symmetric it is a good *toy* example for those 3D functions like shading or slope projections.\n", - "* [Random](#Random---Normal) : \"white-noise\" to be used in `geo.mean` function\n", - " * The visual effect of the mean filter is to *blur* the image. Another important property of the mean filter is that is conserves the total, integrate intensity of the image. Both aspects are clearly visible from a uniform (random) distribution.\n", - " \n", - "#### Real\n", - "* [Austin-DEM](#Austin---DEM) : image providing elevation for the region surrounding the city of Austin, TX, USA; we will use the Austin-DEM data to produce a shaded (`hillshade`) version of the image.\n", - "* [Landsat](#Landsat) : Landsat is a Earth observation mission -- currently in its 8th phase -- providing satellite images in different wavelengths. In particular, we will use two images of the same region but different waveband -- red and near-infrared -- to estimate the vegetation coverage from the `ndvi` method." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Aspect\n", - "* Burrough, P. A., and McDonell, R. A., 1998. *Principles of Geographical Information Systems* (Oxford University Press, New York), pp 406\n", - "\n", - "[Aspect](https://en.wikipedia.org/wiki/Aspect_(geography)) is the orientation of slope, measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "display(Markdown(aspect.__doc__))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# We will use a gaussian surface to help us understand the result of the *aspect* transform.\n", - "\n", - "xarray_gaussian = helper_do_gaussian()\n", - "\n", - "da_aspect = aspect(xarray_gaussian)\n", - "\n", - "hv.Image(da_aspect.data, label='Aspect of Gaussian centered at [0,0]').options(colorbar=True).options(cmap=COLORMAP)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Hillshade\n", - "* Burrough, P. A., and McDonell, R. A., 1998. *Principles of Geographical Information Systems* (Oxford University Press, New York), pp 406\n", - "\n", - "[Hillshade](https://en.wikipedia.org/wiki/Terrain_cartography) is a technique used to visualize terrain as shaded relief, illuminating it with a hypothetical light source. The illumination value for each cell is determined by its orientation to the light source, which is based on slope and aspect." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "display(Markdown(hillshade.__doc__))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Multi-Directional Oblique (MDOW) Hillshade" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "hillshade_agg = hillshade(xarray_gaussian, how='mdow', out_type='data')\n", - "img_shade = hv.Image(hillshade_agg.data, label='Shade (MDOW) from Gaussian')\n", - "img_shade.options(colorbar=True).options(cmap='gray') + img_gauss.options(cmap=COLORMAP)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## NDVI\n", - "* Burrough, P. A., and McDonell, R. A., 1998. *Principles of Geographical Information Systems* (Oxford University Press, New York), pp 406\n", - "\n", - "The Normalized Difference Vegetation Index (NDVI) quantifies vegetation by measuring the difference between near-infrared (which vegetation strongly reflects) and red light (which vegetation absorbs).\n", - "NDVI always ranges from -1 to +1. But there isn’t a distinct boundary for each type of land cover.\n", - "For example, when you have negative values, it’s highly likely that it’s water. On the other hand, if you have a NDVI value close to +1, there’s a high possibility that it’s dense green leaves.\n", - "But when NDVI is close to zero, there isn’t green leaves and it could even be an urbanized area." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "display(Markdown(ndvi.__doc__))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import datashader.transfer_functions as tf\n", - "\n", - "NIR_data = helper_do_gaussian()\n", - "RED_data = helper_do_gaussian_inverse()\n", - "\n", - "da_ndvi = ndvi(NIR_data, RED_data)\n", - "hv.Image(da_ndvi.data).options(colorbar=True).options(cmap='PRGn')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Real data\n", - "\n", - "Here we are going to see how those methods transform real images.\n", - "The images represent the region around Austin, TX in USA, which is mountainous region and make a good example for our `hillshade` method.\n", - "The other data set is composed by two images from the Landsat satellite, the *red* and *near-infrared* filters, which we will use for checking the `ndvi` processing effects." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Austin terrain" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "from datashader import transfer_functions as tf\n", - " \n", - "import rasterio\n", - " \n", - "# Austin-DEM\n", - "austin_dem = '../assets/images/austin_dem.tif'\n", - "\n", - "elevation_data = rasterio.open(austin_dem)\n", - "elevation_xarray = xr.DataArray(elevation_data.read().astype(np.uint16))\n", - "\n", - "xmin = elevation_data.bounds.left\n", - "ymin = elevation_data.bounds.bottom\n", - "xmax = elevation_data.bounds.right\n", - "ymax = elevation_data.bounds.top\n", - "\n", - "w = np.ceil(elevation_data.width / 8.0)\n", - "h = np.ceil(elevation_data.height / 8.0)\n", - "\n", - "cvs = ds.Canvas(plot_width=w, plot_height=h, \n", - " x_range=(xmin,xmax), y_range=(ymin,ymax))\n", - "agg = cvs.raster(elevation_xarray[0])\n", - "hillshade(agg, out_type='data')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "da_elevation = xr.open_rasterio(austin_dem)\n", - "\n", - "# Define Datashader Canvas\n", - "from datashader import Canvas\n", - "from math import ceil\n", - "\n", - "# xmin = elevation_data.bounds.left\n", - "xmin = da_elevation.x.min()\n", - "\n", - "# ymin = elevation_data.bounds.bottom\n", - "ymin = da_elevation.y.min()\n", - "\n", - "# xmax = elevation_data.bounds.right\n", - "xmax = da_elevation.x.max()\n", - "\n", - "# ymax = elevation_data.bounds.top\n", - "ymax = da_elevation.y.max()\n", - "\n", - "# w = np.ceil(elevation_data.width / 8.0)\n", - "w = ceil(len(da_elevation.x) / 8.0)\n", - "\n", - "# h = np.ceil(elevation_data.height / 8.0)\n", - "h = ceil(len(da_elevation.y) / 8.0)\n", - "\n", - "cvs = Canvas(plot_width=w, plot_height=h, x_range=(xmin,xmax), y_range=(ymin,ymax))\n", - "\n", - "xarr = da_elevation.load()[0]\n", - "\n", - "# Rasterize our data according to the Canvas resolution/definitions\n", - "agg = cvs.raster(xarr)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(xarr.data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "agg.data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cvs.raster?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(agg.data)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.hist(elevation_xarray[0].data.flat, bins=100);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Hillshade" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "austin_shade = hillshade(xr.DataArray(data_austin), out_type='data')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from datashader import transfer_functions as tf\n", - "img = tf.shade(austin_shade, cmap='gray', how='linear')\n", - "img" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "hv.Image(austin_shade.data, label='Austin Mountains').options(colorbar=True).options(cmap='gray')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "austin_shade.stack" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Visualize the shade\n", - "cmap = plt.cm.terrain\n", - "rgba = plt.Normalize()(data_austin)\n", - "rgba = cmap(rgba)\n", - "rgba[..., -1] = 1 - 0.5*normalize(austin_shade.data)\n", - "plt.imshow(rgba);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Landsat Red/NIR" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Open the elevetaion terrain file\n", - "#\n", - "import rasterio\n", - "\n", - "# Landsat\n", - "landsat_red = 'assets/images/MERCATOR_LC80210392016114LGN00_B4.TIF'\n", - "landsat_nir = 'assets/images/MERCATOR_LC80210392016114LGN00_B5.TIF'\n", - "\n", - "with rasterio.open(landsat_red) as fp:\n", - " data_landsat_red = fp.read()\n", - "\n", - "with rasterio.open(landsat_nir) as fp:\n", - " data_landsat_nir = fp.read()\n", - "\n", - "data_landsat_red = data_landsat_red[0]\n", - "data_landsat_nir = data_landsat_nir[0]\n", - "print('Image shapes:', data_landsat_red.shape, data_landsat_nir.shape)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "img_red = hv.Image(data_landsat_red, label='Landsat Red filter').options(colorbar=True).options(cmap='Greens')\n", - "img_nir = hv.Image(data_landsat_nir, label='Landsat NIR filter').options(colorbar=True).options(cmap='Blues')\n", - "img_red + img_nir" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Let's cut a region in the middle to have a closer view of the scene\n", - "#\n", - "data_landsat_red = data_landsat_red[2000:5000, 2000:5000]\n", - "data_landsat_nir = data_landsat_nir[2000:5000, 2000:5000]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "img_red = hv.Image(data_landsat_red, label='Landsat Red filter').options(colorbar=True).options(cmap='Greens')\n", - "img_nir = hv.Image(data_landsat_nir, label='Landsat NIR filter').options(colorbar=True).options(cmap='Blues')\n", - "img_red + img_nir" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As we may notice, the images intensities span quite different range of values.\n", - "Let's apply a histogram equalization and normalize the values to have a clearer view of the scene." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from datashader.transfer_functions import eq_hist\n", - "\n", - "data_landsat_red_eq = eq_hist(data_landsat_red)\n", - "data_landsat_nir_eq = eq_hist(data_landsat_nir)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "img_red = hv.Image(data_landsat_red_eq, label='Landsat Red filter').options(colorbar=True).options(cmap='gray')\n", - "img_nir = hv.Image(data_landsat_nir_eq, label='Landsat NIR filter').options(colorbar=True).options(cmap='gray')\n", - "img_red + img_nir" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data_landsat_ndvi = ndvi(xr.DataArray(data_landsat_nir_eq), xr.DataArray(data_landsat_red_eq))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "hv.Image(data_landsat_ndvi.data).options(colorbar=True).options(cmap='PRGn')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The output of *NDVI* ranges from [-1,+1], where `-1` means more \"Red\" radiation while `+1` means more \"NIR\" radiation as observed by Landsat. A higher NIR radiation trace the presence of vegetation." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { - "kernelspec": { - "display_name": "Python 2", - "language": "python", - "name": "python2" - }, "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.15" + "pygments_lexer": "ipython3" } }, "nbformat": 4, diff --git a/tox.ini b/tox.ini index 95989229e..3841f0787 100644 --- a/tox.ini +++ b/tox.ini @@ -63,7 +63,6 @@ norecursedirs = doc .git dist build _build .ipynb_checkpoints # TODO: skipping everything but getting_started/2 for now nbsmoke_skip_run = .*nyc_taxi-nongeo\.ipynb$ .*streaming-aggregation\.ipynb$ - .*geo_transfer_functions.ipynb$ .*topics/.*\.ipynb$ .*user_guide/.*\.ipynb$ .*getting_started/1_Introduction\.ipynb$