Using Geoserver to load data on the map

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

Make sure you have the geoserver VM running

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

Display geoserver status

This should ensure the client can successfully connect to your VM, if you do not see the Geoserver 'Status' page then something is wrong and the rest of the notebook may not function correctly.


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))

Get the data from S3

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

Adding an RGB layer to the map

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]:
<promise.promise.Promise at 0x7fb5c7855f10>

In [3]:
from geonotebook.wrappers import RasterData

rd = RasterData('data/L57.Globe.month09.2010.hh09vv04.h6v1.doy247to273.NBAR.v3.0.tiff')
rd


Out[3]:
<geonotebook.wrappers.RasterData at 0x7fb5c758b1d0>

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]:
<promise.promise.Promise at 0x7fb5c758bf10>

In [8]:
M.layers.annotation.points[0].data.next()


Out[8]:
(<SimpleLayer('L57.Globe.month09.2010.hh09vv04.h6v1.doy247to273.NBAR.v3.0')>,
 array([-9999., -9999., -9999.], dtype=float32))

In [9]:
from geonotebook.vis.ktile.utils import get_layer_vrt
print get_layer_vrt(M.layers[0])


<VRTDataset rasterXSize="7243" rasterYSize="1900">
    <SRS>EPSG:4326</SRS>
    <GeoTransform>-123.062293136, 0.000751761084359, 0.0, 48.5714292392, 0.0, -0.000751761084359</GeoTransform>
    <VRTRasterBand dataType="Float32" band="1">
        <NoDataValue>-9999.0</NoDataValue>
        <ColorInterp>Red</ColorInterp>
        <ComplexSource>
            <SourceFilename relativeToVRT="0">/home/kotfic/src/geonotebook/notebooks/data/L57.Globe.month09.2010.hh09vv04.h6v1.doy247to273.NBAR.v3.0.tiff</SourceFilename>
            <SourceBand>1</SourceBand>
            <ScaleRatio>2.550000e+02</ScaleRatio>
            <NODATA>-9999.0</NODATA>
        </ComplexSource>
    </VRTRasterBand>
    <VRTRasterBand dataType="Float32" band="2">
        <NoDataValue>-9999.0</NoDataValue>
        <ColorInterp>Green</ColorInterp>
        <ComplexSource>
            <SourceFilename relativeToVRT="0">/home/kotfic/src/geonotebook/notebooks/data/L57.Globe.month09.2010.hh09vv04.h6v1.doy247to273.NBAR.v3.0.tiff</SourceFilename>
            <SourceBand>2</SourceBand>
            <ScaleRatio>2.550000e+02</ScaleRatio>
            <NODATA>-9999.0</NODATA>
        </ComplexSource>
    </VRTRasterBand>
    <VRTRasterBand dataType="Float32" band="3">
        <NoDataValue>-9999.0</NoDataValue>
        <ColorInterp>Blue</ColorInterp>
        <ComplexSource>
            <SourceFilename relativeToVRT="0">/home/kotfic/src/geonotebook/notebooks/data/L57.Globe.month09.2010.hh09vv04.h6v1.doy247to273.NBAR.v3.0.tiff</SourceFilename>
            <SourceBand>3</SourceBand>
            <ScaleRatio>2.550000e+02</ScaleRatio>
            <NODATA>-9999.0</NODATA>
        </ComplexSource>
    </VRTRasterBand>
</VRTDataset>

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])

Adding a single band Layer

Adding a single band layer uses the same M.add_layer(...) interface. Keep in mind that several of the styling options are slightly different. By default single band rasters are rendered with a default mapping of colors to band values.


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]:
<promise.promise.Promise at 0x7fc26ca63950>

What can I do with this data?

We will address the use of annotations for analysis and data comparison in a separate notebook. For now Let's focus on a small agricultural area north of I-90:


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))

NDVI Segmentation analysis

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)