from bokeh.plotting import figure, show, output_notebook
import json
from bokeh import palettes
from bokeh import plotting
import numpy as np
import ogr
import gdal
import pandas as pd
from bokeh.models import HoverTool, TapTool
from bokeh.models import ColumnDataSource
from bokeh.models import Patches
from bokeh.models.widgets import Tabs, Panel
from tools import getAttributes
import pyproj
import osr
inraster = r'D:\Temp\gmted_average_100km_proj.tif'
ds = gdal.Open(inraster)
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
ds = None
cr = osr.SpatialReference()
cr.ImportFromWkt(proj)
inproj = cr.ExportToProj4()
ysize, xsize = data.shape
ulx = gt[0]
uly = gt[3]
lrx = ulx + xsize * gt[1]
lry = uly + ysize * gt[5]
xres = gt[1]
yres = gt[5]
topleft_xs = [gt[0] + gt[1] * i for i in range(xsize)]
topleft_ys = [gt[3] + gt[5] * i for i in range(ysize)]
tl = [(x, y) for x in topleft_xs for y in topleft_ys]
#tl = list(zip(topleft_xs, ))
tr = [(xy[0] + xres, xy[1]) for xy in tl]
br = [(xy[0] + xres, xy[1] + yres) for xy in tl]
bl = [(xy[0], xy[1] + yres) for xy in tl]
polys = list(zip(tl, tr, br, bl, tl))
len(tl)
1620
def toLatLon(coords):
coords = list(zip(*coords))
xs = coords[0]
ys = coords[1]
crobj = pyproj.Proj(inproj)
lons, lats = crobj(xs, ys, inverse=True)
return lons, lats
lonlats = list(map(toLatLon, polys))
len(lonlats)
1620
xs = [xs for xs, ys in lonlats]
ys = [ys for xs, ys in lonlats]
df = pd.DataFrame({'value': data.T.ravel().tolist(),
'xs': xs,
'ys': ys})
color_key = 'value'
cmap = palettes.Blues9[::-1]
vmax = df[color_key].max() / 2.
n = len(cmap)-1
df['color'] = [cmap[min(int(round(val / vmax * n)), n)] for val in df[color_key]]
source = ColumnDataSource(df)
dfcountries = getAttributes(r'D:\Temp\ne_110m_admin_0_countries_lakes_subset.shp', include_geom=True, geom_as_bokeh=True)
cntry_source = ColumnDataSource(dfcountries)
def make_map(source, fill_color='color'):
polys = Patches(xs='xs',
ys='ys',
fill_color=fill_color,
line_color='black',
fill_alpha=.6,
line_width=.15
)
countries = Patches(xs='xs',
ys='ys',
fill_color='#eeeeee',
line_color='black',
fill_alpha=.5,
line_alpha=.5,
line_width=.5
)
xrng = [-140, -55]
yrng = [15, 70]
ratio = (xrng[0] - xrng[1]) / (yrng[0] - yrng[1])
width = 800
height = int(width / ratio)
p = figure(title='Elevation',
plot_width=width,
plot_height=height,
x_range=xrng,
y_range=yrng,
toolbar_location=None,
outline_line_color='#FFFFFF',
min_border=0,
tools='box_zoom, pan, wheel_zoom')
p.add_glyph(cntry_source, countries)
p.add_glyph(source, polys)
p.grid.grid_line_color = None
return p
p = make_map(source)
show(p)