<jupyter_start><jupyter_text>Datashader is a general-purpose tool for rasterizing (and re-rasterizing) data of many different types. To make it easier to apply this general functionality to the particular domain of geoscience, Datashader provides a few geospatial-specific utilities as well:* [Project points](Project-points)* [Generate terrain](Generate-terrain)* [Hillshade](Hillshade)* [Slope](Slope)* [Aspect](Aspecct)* [Bump map](Bump-map)* [NDVI](NDVI)* [Mean](Mean)This notebook explains each of these topics in turn. See also [GeoViews](http://geoviews.org), which is designed to work with Datashader to provide a large range of additional geospatial functionality.  Project pointsYou can use [GeoViews](http://geoviews.org) or the underlying [pyproj/proj.4](https://jswhit.github.io/pyproj) libraries to perform arbitrary projections to and from a large number of different coordinate reference systems. However, for the common case of wanting to view data with latitude and longitude coordinates on top of a Web Mercator tile source such as Google Maps or OpenStreetMap, Datashader provides a simple self-contained utility `lnglat_to_meters(longitude, latitude)` to project your data once, before visualization.  For instance, if you have a dataframe with some latitude and longitude points stretching from San Diego, California to Bangor, Maine:<jupyter_code>import numpy as np, pandas as pd
from datashader.utils import lnglat_to_meters

San_Diego = 32.715, -117.1625
Bangor = 44.8, -68.8
n = 20

df = pd.DataFrame(dict(longitude = np.linspace(San_Diego[1], Bangor[1], n),
                       latitude  = np.linspace(San_Diego[0], Bangor[0], n)))<jupyter_output><empty_output><jupyter_text>Then you can create new columns (or overwrite old ones) with the projected points in meters from the origin (Web Mercator coordinates):<jupyter_code>df.loc[:, 'x'], df.loc[:, 'y'] = lnglat_to_meters(df.longitude,df.latitude)
df.tail()<jupyter_output><empty_output><jupyter_text>The new x and y coordinates aren't very useful for humans to read, but they can now be overlaid directly onto web map sources, which are labeled with latitude and longitude appropriately by Bokeh but are actually in Web Mercator coordinates internally:<jupyter_code>import holoviews as hv, geoviews.tile_sources as gts
from holoviews.operation.datashader import datashade, spread
hv.extension('bokeh')

gts.EsriImagery() * spread(datashade(hv.Points(df, ['x', 'y']), cmap="white"), px=3)<jupyter_output><empty_output><jupyter_text>If you are using GeoViews, you can get the same effect by calling [gv.operation.project](http://geoviews.org/user_guide/Projections.htmlExplicitly-projecting).  With GeoViews, you can also declare your object to be in lon,lat coordinates natively (`from cartopy import crs ; gv.Points(df, ['longitude', 'latitude'], crs=crs.PlateCarree())`) and let GeoViews then reproject the points as needed, but dynamic reprojection will be much slower for interactive use than projecting them in bulk ahead of time.  Generate Terrain DataThe rest of the geo-related functions focus on raster data (or rasterized data, after a previous Datashader step that returns an Xarray object). To demonstrate using these raster-based functions, let's generate some fake terrain as an elevation raster:<jupyter_code>import numpy as np, datashader as ds, datashader.geo as dsgeo
from datashader.transfer_functions import shade, stack
from datashader.colors import Elevation

W = 1000
H = 750

canvas = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6))
terrain = dsgeo.generate_terrain(canvas)

shade(terrain, cmap=['black', 'white'], how='linear')<jupyter_output><empty_output><jupyter_text>The grayscale value above shows the elevation linearly in intensity (with the large black areas indicating low elevation), but it will look more like a landscape if we map the lowest values to colors representing water, and the highest to colors representing mountaintops:<jupyter_code>shade(terrain, cmap=Elevation, how='linear')<jupyter_output><empty_output><jupyter_text>Hillshade[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.<jupyter_code>illuminated = dsgeo.hillshade(terrain)

shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')<jupyter_output><empty_output><jupyter_text>You can combine hillshading with elevation colormapping to convey differences in terrain with elevation:<jupyter_code>stack(shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear'),
      shade(terrain,     cmap=Elevation,         alpha=128, how='linear'))<jupyter_output><empty_output><jupyter_text>Slope[Slope](https://en.wikipedia.org/wiki/Slope) is the inclination of a surface. In geography, *slope* is amount of change in elevation of a terrain regarding its surroundings.Datashader's slope function returns slope in degrees.  Below we highlight areas at risk for avalanche by looking at [slopes around 38 degrees](http://wenatcheeoutdoors.org/2016/04/07/avalanche-abcs-for-snowshoers/).<jupyter_code>risky = dsgeo.slope(terrain)
risky.data = np.where(np.logical_and(risky.data > 25, risky.data < 50), 1, np.nan)

stack(shade(terrain,     cmap=['black', 'white'], how='linear'),
      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),
      shade(risky,       cmap='red',              how='linear', alpha=200))<jupyter_output><empty_output><jupyter_text>Aspect[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.Below, we look to find slopes that face close to North.<jupyter_code>north_faces = dsgeo.aspect(terrain)
north_faces.data = np.where(np.logical_or(north_faces.data > 350 ,
                                          north_faces.data < 10), 1, np.nan)

stack(shade(terrain,     cmap=['black', 'white'], how='linear'),
      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),
      shade(north_faces, cmap=['aqua'],           how='linear', alpha=100))<jupyter_output><empty_output><jupyter_text>NDVIThe Normalized Difference Vegetation Index (NDVI) quantifies vegetation by measuring the difference between near-infrared (which vegetation strongly reflects) and red light (which vegetation absorbs).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.But when NDVI is close to zero, there isn’t green leaves and it could even be an urbanized area. The output of *NDVI* ranges from [-1,+1], where `-1` means more "Red" radiation while `+1` means more "NIR" radiation.Below, we simulate the red and near-infrared bands using `datashader.perlin` random noise with different seeds and frequencies.  Green areas should be those > 0, where higher NDVI values would indicate vegetation.<jupyter_code>near_infrared_band = dsgeo.perlin(W, H, freq=(6,   6), seed=1)
red_band           = dsgeo.perlin(W, H, freq=(10, 10), seed=2)
vegetation_index   = dsgeo.ndvi(near_infrared_band, red_band)

shade(vegetation_index, cmap=['purple','black','green'], how='linear')<jupyter_output><empty_output><jupyter_text>BumpBump mapping is a cartographic technique that can be used to create the appearance of trees or other land features.`dsgeo.bump` will produce a bump aggregate for adding detail to the terrain.  In this case, we will pretend the bumps are trees, and shade them with green.<jupyter_code>from functools import partial

def heights(locations, min_val, max_val, height):
    out = np.zeros(len(locations))
    for i, (x, y) in enumerate(locations):
        val = terrain.data[y, x]
        out[i] = height if (val > min_val and val < max_val) else 0
    return out

T = 200000

trees  = dsgeo.bump(W, H, count=T//3, height_func=partial(heights, min_val=50,   max_val=500,  height=10))
trees += dsgeo.bump(W, H, count=T,    height_func=partial(heights, min_val=500,  max_val=2000, height=20))
trees += dsgeo.bump(W, H, count=T//3, height_func=partial(heights, min_val=2000, max_val=3000, height=10))

tree_colorize = trees.copy()
tree_colorize.data[tree_colorize.data == 0] = np.nan

stack(shade(terrain + trees,                  cmap=['black', 'white'], how='linear'),
      shade(dsgeo.hillshade(terrain + trees), cmap=['black', 'white'], how='linear', alpha=128),
      shade(tree_colorize,                    cmap='limegreen',        how='linear'))<jupyter_output><empty_output><jupyter_text>MeanThe `datashader.mean` function will smooth a given aggregate by using a 3x3 mean convolution filter. Optional parameters include `passes`, which is used to run the mean filter multiple times, and also `excludes` which are values that will not be modified by the mean filter.Just for fun, let's add a coastal vignette to give out terrain scene a bit more character. Notice the water below now has a nice coastal gradient which adds some realism to our scene.<jupyter_code>LAND_CONSTANT = 50.0

water = terrain.copy()
water.data = np.where(water.data > 0, LAND_CONSTANT, 0)
water = dsgeo.mean(water, passes=10, excludes=[LAND_CONSTANT])
water.data[water.data == LAND_CONSTANT] = np.nan

stack(shade(terrain,                  cmap=['black', 'white'], how='linear'),
      shade(water,                    cmap=['aqua',  'white']),
      shade(dsgeo.hillshade(terrain), cmap=['black', 'white'], how='linear', alpha=128))<jupyter_output><empty_output><jupyter_text>ConclusionWe've now seen a bunch of datashader's `geo` helper functions for working with elevation rasters.Let's make our final archipelago scene by stacking `terrain`, `water`, `hillshade`, and `tree_highlights` together into one output image:<jupyter_code>stack(shade(terrain + trees,                  cmap=Elevation,          how='linear'),
      shade(water,                            cmap=['aqua','white']),
      shade(dsgeo.hillshade(terrain + trees), cmap=['black', 'white'], how='linear', alpha=128),
      shade(tree_colorize,                    cmap='limegreen',        how='linear'))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1069816575434553, 'fraction_numerical': 0.01613367905502737, 'mean_word_length': 5.1294879340788695, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 58, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '635466', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3393, 'n_tokens_neox': 3141, 'n_words': 1054, 'path': 'examples/user_guide/8_Geography.ipynb', 'repo_name': 'seanlane/datashader', 'size': '10413'}