ORTHOGONAL POLYNOMIAL DENSITY ESTIMATION

In [1]:

    

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.colors import ListedColormap
from mpl_toolkits.axes_grid1 import make_axes_locatable
from uuid import uuid4

from lpde.geometry import WidthOf, Window, PointAt, BoundingBox, Mapper, Grid
from lpde.estimators import ParallelEstimator
from lpde.estimators.datatypes import Event, Degree, Action
from lpde.producers import MockParams
from lpde.producers.distributions import gaussian
from lpde.visualizers import Visualize


    

In [2]:

    

%matplotlib qt5


    

In [3]:

    

legendre_width = WidthOf(1.8)

center = PointAt(51.375, 35.675)
window = Window(0.55, 0.35)
bounds = BoundingBox(center, window)

mapper = Mapper(bounds, legendre_width)

degree = Degree(20, 20)
params = MockParams(20, 100, gaussian)
demand = ParallelEstimator(degree, mapper, params)

action = Action.ADD
point = PointAt(51.2, 35.6)
event = Event(uuid4(), action, point)

grid = Grid(100, 100)


    

In [4]:

    

demand.controller.start(1, 1.0)


    

    




/home/georg/Documents/Python/LPDE/lpde/estimators/parallel/minimizer.py:117: RuntimeWarning: overflow encountered in double_scalars
  return self.__neg_log_l(c[1:]) + c[0]*self.__norm(c[1:])
/home/georg/Documents/Python/LPDE/lpde/estimators/parallel/minimizer.py:125: RuntimeWarning: overflow encountered in square
  return -log(square(c.dot(self.__phi_ijn))).sum()
/home/georg/Documents/Python/LPDE/lpde/estimators/parallel/minimizer.py:121: RuntimeWarning: overflow encountered in multiply
  self.__grad_c[1:] = self.__grad_neg_log_l(c[1:]) + 2.0*c[0]*c[1:]

In [5]:

    

demand.controller.alive


    

    
Out[5]:





{'datagate': True, 'minimizers': (True,), 'producer': True, 'smoother': True}

In [6]:

    

demand.controller.open


    

    
Out[6]:





{'coefficients': True, 'events': True, 'points': True}

In [7]:

    

demand.controller.qsize


    

    
Out[7]:





{'coefficients': 0, 'events': 0, 'points': 26}

In [ ]:

    

demand.update_with(event)


    

In [11]:

    

demand.controller.stop()


    

In [ ]:

    

visulization = Visualize(demand)


    

In [ ]:

    

animation = visulization.show(cartopy=True, zoom=11)


    

In [8]:

    

import cartopy.crs as cpcrs
import cartopy.io.img_tiles as cpimgt


    

In [9]:

    

mpl_cmap = plt.cm.viridis
new_cmap = mpl_cmap(np.arange(mpl_cmap.N))
new_cmap[:, -1] = np.linspace(0, 1, mpl_cmap.N)
new_cmap = ListedColormap(new_cmap)


    

In [29]:

    

fig = plt.figure()

osm = cpimgt.OSM()

iax = plt.axes(projection=osm.crs)
iax.set_extent(bounds.extent)
iax.add_image(osm, 10)
iax.text(-0.05, 0.50, 'latitude', rotation='vertical', va='center', transform=iax.transAxes)
iax.text(0.5, -0.05, 'longitude', rotation='horizontal', ha='center', transform=iax.transAxes)

contour = iax.imshow(demand.on_grid,
                     transform=cpcrs.PlateCarree(),
                     cmap=new_cmap,
                     extent=bounds.extent,
                     origin='lower',
                     interpolation='bilinear',
                     animated=True)

def initialize_figure():
    return contour,

def update_figure(*args):
    data = demand.on_grid
    contour.set_data(data)
    contour.set_clim(data.min(), data.max())
    return contour,

animation = FuncAnimation(fig, update_figure, init_func=initialize_figure, interval=50, blit=True)

plt.show()


    

In [ ]:

    

fig = plt.figure()

osm = cpimgt.OSM()

iax = plt.axes(projection=osm.crs)
iax.set_extent(bounds.extent)
iax.add_image(osm, 10)
iax.text(-0.05, 0.50, 'latitude', rotation='vertical', va='center', transform=iax.transAxes)
iax.text(0.5, -0.05, 'longitude', rotation='horizontal', ha='center', transform=iax.transAxes)

contour = iax.imshow(demand.on_grid,
                     transform=cpcrs.PlateCarree(),
                     cmap=new_cmap,
                     extent=bounds.extent,
                     origin='lower',
                     interpolation='bilinear',
                     animated=True)

def initialize_figure():
    return contour,

def update_figure(*args):
    data = demand.on_grid
    contour.set_data(data)
    contour.set_clim(data.min(), data.max())
    return contour,

animation = FuncAnimation(fig, update_figure, init_func=initialize_figure, interval=50, blit=True)

plt.show()


    

In [ ]:

    

fig_size_x = 8

fig = plt.figure(figsize=(fig_size_x, fig_size_x*bounds.aspect))
ax = plt.axes()
axlabel = ax.set(xlabel='longitude', ylabel='lattitude')

contour = ax.imshow(demand.on_grid,
                    cmap='viridis',
                    extent=bounds.x_range+bounds.y_range,
                    origin='lower',
                    animated=True)

divider = make_axes_locatable(ax)
cbax = divider.append_axes('right', size='5%', pad=0.1)
cbar = plt.colorbar(contour, cax=cbax, label='demand')

def initialize_figure():
    return contour,

def update_figure(*args):
    data = demand.on_grid
    contour.set_data(data)
    contour.set_clim(data.min(), data.max())
    return contour,

animation = FuncAnimation(fig, update_figure, interval=100, blit=False)

fig.tight_layout()
fig.show()


    

In [ ]:

    

bokeh_extent = iax.get_extent()
bokeh_x_range = bokeh_extent[:2]
bokeh_y_range = bokeh_extent[2:]
bokeh_x = bokeh_extent[0]
bokeh_y = bokeh_extent[2]
bokeh_dx = bokeh_extent[1] - bokeh_extent[0]
bokeh_dy = bokeh_extent[3] - bokeh_extent[2]

image = contour.make_image('Agg')[0]


    

In [ ]:

    

from bokeh.io import output_file, show
from bokeh.plotting import figure
from bokeh.tile_providers import CARTODBPOSITRON
from bokeh.palettes import Viridis256


fig = figure(tools='pan, wheel_zoom', x_range=bokeh_x_range, y_range=bokeh_y_range)
fig.axis.visible = False
fig.add_tile(CARTODBPOSITRON)
fig.image_rgba([image], x=bokeh_x, y=bokeh_y, dw=bokeh_dx, dh=bokeh_dy)
output_file("stamen_toner_plot.html")
show(fig)


    

In [ ]: