In this notebook we'll take a look at using Geoserver to render raster data to the map. Geoserver is an open source server for sharing geospatial data. It includes a tiling server which the GeoJS map uses to render data efficiently to the map for visualization. Geonotebook comes with a vagrant virtual machine for hosting a local instance of Geoserver. This instance can be used for testing geonotebook. To use it simply install vagrant using your system package manager, in a checked out copy of the source code go to the devops/geoserver/
folder and run vagrant up
In [1]:
%matplotlib inline
from matplotlib import pylab as plt
The following cell will check whether or not your have a running instance of the geoserver virtual machine available. The following cell should show text to the effect of:
Current machine states:
geoserver running (virtualbox)
The VM is running. To stop this VM, you can run `vagrant halt` to
shut it down forcefully, or you can run `vagrant suspend` to simply
suspend the virtual machine. In either case, to restart it again,
simply run `vagrant up`.
If it does not show the geoserver machine in a state of running
You can load the machine by going to ../devops/geoserver/
and running vagrant up
In [ ]:
!cd ../devops/geoserver && vagrant status
In [ ]:
from IPython.core.display import display, HTML
from geonotebook.config import Config
geoserver = Config().vis_server
display(HTML(geoserver.c.get("/about/status").text))
Next get some sample data from S3. This GeoTiff represents NBAR data for September from 2010 covering a section of Washington states Glacier National Park. It is aproximately 200Mb and may take some time to download from Amazon's S3.
The tiff itself has been slightly transformed from its original HDF dataset. In particular it only has 4 bands (R,G,B & NDVI) and includes some geotiff tags with band statistics.
In [ ]:
!curl -o /tmp/L57.Globe.month09.2010.hh09vv04.h6v1.doy247to273.NBAR.v3.0.tiff http://golden-tile-geotiffs.s3.amazonaws.com/L57.Globe.month09.2010.hh09vv04.h6v1.doy247to273.NBAR.v3.0.tiff
Here we add our first data layer to the map. To do this we use a RasterData
object imported from the geonotebook.wrappers package. By default RasterData objects read tiffs using the rasterio
library. RasterData objects are designed to provide a consistent API to raster data across a number of different readers and systems. We will use the add_layer
function to add the RasterData object to the map.
In [2]:
# Set the center of the map to the location the data
M.set_center(-120.32, 47.84, 7)
Out[2]:
In [3]:
from geonotebook.wrappers import RasterData
rd = RasterData('data/L57.Globe.month09.2010.hh09vv04.h6v1.doy247to273.NBAR.v3.0.tiff')
rd
Out[3]:
To add the layer we call M.add_layer
passing in a subset of the raster data set's bands. In this case we index into rd
with the list [1, 2, 3]. This actually returns a new RasterData object with only three bands available (in this case bands 1, 2 and 3 corrispond to Red, Green and Blue). When adding layers you can only add a layer with either 3 bands (R,G,B) or one band (we'll see a one band example in a moment).
In [6]:
M.add_layer(rd[1, 2, 3], opacity=1.0)
Out[6]:
In [8]:
M.layers.annotation.points[0].data.next()
Out[8]:
In [9]:
from geonotebook.vis.ktile.utils import get_layer_vrt
print get_layer_vrt(M.layers[0])
This should have added an RGB dataset to the map for visualization. You can also see what layers are available via the M.layers
attribute.
In [ ]:
M.layers
The dataset may appear alarmingly dark. This is because the data itself is not well formated. We can see this by looking at band min and max values:
In [ ]:
print("Color Min Max")
print("Red: {}, {}".format(rd[1].min, rd[1].max))
print("Green: {}, {}".format(rd[2].min, rd[2].max))
print("Blue: {}, {}".format(rd[3].min, rd[3].max))
R,G,B values should be between 0 and 1. We can remedy this by changing some of the styling options that are available on the layers including setting an interval for scaling our data, and setting a gamma to brighten the image.
First we'll demonstrate removing the layer:
In [ ]:
M.remove_layer(M.layers[0])
Then we can re-add the layer with a color interval of 0 to 1.
In [ ]:
M.add_layer(rd[1, 2, 3], interval=(0,1))
We can also brighten this up by changing the gamma.
Note We don't have to remove the layer before updating it's options. Calling M.add_layer(...)
with the same rd
object will simply replace any existing layer with the same name. By default the layer's name is inferred from the filename.
In [ ]:
M.add_layer(rd[1, 2, 3], interval=(0,1), gamma=0.5)
Finally, let's add a little opacity to layer so we can see some of the underlying base map features.
In [ ]:
M.add_layer(rd[1, 2, 3], interval=(0,1), gamma=0.5, opacity=0.75)
In [ ]:
# Remove the layer before moving on to the next section
M.remove_layer(M.layers[0])
In [ ]:
M.add_layer(rd[4])
You may find this colormap a little aggressive, in which case you can replace the colormap with any of the built in matplotlib colormaps:
In [ ]:
cmap = plt.get_cmap('winter', 10)
M.add_layer(rd[4], colormap=cmap, opacity=0.8)
Including custom color maps as in this example. Here we create a linear segmented colormap that transitions from Blue to Beige to Green. When mapped to our NDVI band data -1 will appear blue, 0 will appear beige and 1 will appear green.
In [6]:
from matplotlib.colors import LinearSegmentedColormap
# Divergent Blue to Beige to Green colormap
cmap =LinearSegmentedColormap.from_list(
'ndvi', ['blue', 'beige', 'green'], 20)
# Add layer with custom colormap
M.add_layer(rd[4], colormap=cmap, opacity=0.8, min=-1.0, max=1.0)
Out[6]:
In [ ]:
M.set_center(-119.25618502500376, 47.349300631765104, 11)
Go ahead and start a rectangular annotation (Second button to the right of the 'CellToolbar' button - with the square icon).
Please annotate a small region of the fields.
We can access this data from from the annotation's data
attribute. We'll cover exactly what is going on here in another notebook.
In [ ]:
layer, data = next(M.layers.annotation.rectangles[0].data)
data
As a sanity check we can prove the data is the region we've selected by plotting the data with matplotlib's imshow function:
Note The scale of the matplotlib image may seem slightly different than the rectangle you've selected on the map. This is because the map is displaying in Web Mercator projection (EPSG:3857) while imshow is simply displaying the raw data, selected out of the geotiff (you can think of it as being in a 'row', 'column' projection).
In [ ]:
import numpy as np
fig, ax = plt.subplots(figsize=(16, 16))
ax.imshow(data, interpolation='none', cmap=cmap, clim=(-1.0, 1.0))
Once we have this data we can run arbitrary analyses on it. In the next cell we use a sobel filter and a watershed transformation to generate a binary mask of the data. We then use an implementation of marching cubes to vectorize the data, effectively segmenting green areas (e.g. fields) from surrounding areas.
This next cell requires both scipy and scikit-image. Check your operating system documentation for how best to install these packages.
In [ ]:
# Adapted from the scikit-image segmentation tutorial
# See: http://scikit-image.org/docs/dev/user_guide/tutorial_segmentation.html
import numpy as np
from skimage import measure
from skimage.filters import sobel
from skimage.morphology import watershed
from scipy import ndimage as ndi
THRESHOLD = 20
WATER_MIN = 0.2
WATER_MAX = 0.6
fig, ax = plt.subplots(figsize=(16, 16))
edges = sobel(data)
markers = np.zeros_like(data)
markers[data > WATER_MIN] = 2
markers[data > WATER_MAX] = 1
mask = (watershed(edges, markers) - 1).astype(bool)
seg = np.zeros_like(mask, dtype=int)
seg[~mask] = 1
# Fill holes
seg = ndi.binary_fill_holes(seg)
# Ignore entities smaller than THRESHOLD
label_objects, _ = ndi.label(seg)
sizes = np.bincount(label_objects.ravel())
mask_sizes = sizes > THRESHOLD
mask_sizes[0] = 0
clean_segs = mask_sizes[label_objects]
# Find contours of the segmented data
contours = measure.find_contours(clean_segs, 0)
ax.imshow(data, interpolation='none', cmap=cmap, clim=(-1.0, 1.0))
ax.axis('tight')
for n, contour in enumerate(contours):
ax.plot(contour[:, 1], contour[:, 0], linewidth=4)