4/21/201787.8k
adeaner
Public

Hackathon April 2017 Muddy Waters V2

Raster mask of water by sampling OSM water features in vector services and deferred processing and tile serving from PDAP. This RADD Notebook is an example client showing the results.


DigitalGlobe (Required Remote Kernel)
Some of the data for this notebook was provided by the DigitalGlobe remote kernel. You will need access to this remote kernel to be able to make full use of this notebook.

# %matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from gbdx_auth import gbdx_auth
import math
from rasterio.io import MemoryFile
from mpl_toolkits.axes_grid1 import make_axes_locatable
import PIL

from IPython.display import display
import matplotlib.image as mpimg


session = gbdx_auth.get_session()
  
pdap_graph_url =  'http://idahoapitest.geobigdata.io/v1/graph'
    
# TODO: send request to muddy waters, get graph id

id = "57134443-dc3b-4b20-af57-f7706130df3d"

pdap_metadata_url =  'http://idahoapitest.geobigdata.io/v1/metadata/idaho-virtual/'

image = session.get(pdap_metadata_url + id + '/Invert/image.json').json()
georef = session.get(pdap_metadata_url + id + '/Invert/georeferencing.json').json()

image = session.get(pdap_metadata_url + id + '/Invert/image.json').json()
tileXMin = image['tileXOffset']
tileXMax = image['tileXOffset'] + min(image['numXTiles'], 4) - 1
tileYMin = image['tileYOffset']
tileYMax = image['tileYOffset'] + min(image['numYTiles'], 4) - 1

print("%d, %d" % (tileXMin, tileXMax))
print("%d, %d" % (tileYMin, tileYMax))
0, 3
34, 37
# Take red, green blue channels from a multispectral image and convert to uint8 RGB
def display_raw(im, band=0):
    nbands, x, y = im.shape
    if band > nbands:
        raise IndexError("Specified bound outside of image bound range")
    fig, ax = plt.subplots()
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='5%', pad=0.05)
    imx = ax.imshow(im[band,:,:].squeeze(), cmap='bone')
    fig.colorbar(imx, cax=cax, orientation='vertical')
    ax.set_axis_off()
    ax.grid(False)
    return ax
    
pngs = []
pdap_png_url = 'http://idahoapitest.geobigdata.io/v1/tile/idaho-virtual/'
for tileY in range(tileYMin, tileYMax):
    for tileX in range(tileXMin, tileXMax):
        print("%d, %d" % (tileX, tileY))
        png_resp = session.get(pdap_png_url + id + '/Invert/' + str(tileX) + '/' + str(tileY) + '.png')
        with MemoryFile(png_resp.content) as memfile:
            with memfile.open() as dataset:
                data = dataset.read()
        pngs.append(data)
        
imgs = []
for tileY in range(tileYMin, tileYMax):
    for tileX in range(tileXMin, tileXMax):
        tileIdx = (tileY - tileYMin) * (tileXMax - tileXMin) + (tileX - tileXMin)
        print(tileIdx)
        print pngs[tileIdx].shape
        img = display_raw(pngs[tileIdx])
        imgs.append(img)
        
def display_all():
    for img in imgs:
        ''''''''
        
        display(imgs)
0, 34
1, 34
2, 34
0, 35
1, 35
2, 35
0, 36
1, 36
2, 36
0
(1, 256, 256)
1
(1, 256, 256)
2
(1, 256, 256)
3
(1, 256, 256)
4
(1, 256, 256)
5
(1, 256, 256)
6
(1, 256, 256)
7
(1, 256, 256)
8
(1, 256, 256)