source: src/django_gheat/website/tile.py@ 9172

#!/usr/bin/env python
#
# Hack to show image generation realtime, sample tile server implementation.
#
# Rick van der Zwet <info@rickvanderzwet.nl>
from django.core.management import setup_environ
from django.db.models import Max
from django.http import HttpResponse
from gheat.models import *
import logging
import pygame
import sys
import tempfile

# Rending with PIL and computation with numpy has proven to be to slow to be
# usable, but is still in here for refence purposes.
try:
  from PIL import Image
  import ImageDraw
  import numpy as np
except ImportError:
  pass

logging.basicConfig(level=logging.WARNING)
log = logging.getLogger('tile')

class PyGamePicture():
  """ Basic PyGame class, allowing simple image manipulations """
  def __init__(self, method, size):
    self.surf = pygame.Surface(size,flags=pygame.SRCALPHA)

  def center_crop(self,size):
    """ Resize to make centered rectange from image """
    new_surf = pygame.Surface(size, flags=pygame.SRCALPHA)
    curr_size = self.surf.get_size()
    new_surf.blit(self.surf,(0,0),
      ((curr_size[0] - size[0]) / 2, (curr_size[1] - size[1]) / 2, size[0], size[1]))
    self.surf = new_surf

  def write(self, fh,format='png'):
    # XXX: How to get a PNG stream directly to the output
    f = tempfile.NamedTemporaryFile(suffix=format)
    pygame.image.save(self.surf,f.name)
    f.seek(0)
    fh.write(f.read())


  def add_circle(self, center, radius, colour=(255,0,0), transparancy=0):
    """ 
    Hack to add lineair gradient circles and merge with the parent. The
    transparancy can be configured to make the circles to fade out in the
    beginning
    """
    # Make calculations and ranges a whole bunch more easy
    radius = int(math.ceil(radius))

    new_surf = pygame.Surface(self.surf.get_size(),flags=pygame.SRCALPHA)
    alpha_per_radius = float(2.55 * (100 - transparancy)) / radius
    for r in range(radius,1,-1):
      alpha = int((radius - r) * alpha_per_radius)
      pygame.draw.circle(new_surf,colour + (alpha,),center,r,0)
    self.surf.blit(new_surf,(0,0),special_flags=pygame.BLEND_RGBA_MAX)


class PILPicture():
  """ Basic PIL class, allowing simple image manipulations """
  im = None
  def __init__(self, method, size):
    self.im = Image.new(method, size)
    self.data = np.array(self.im)

  def write(self,fh,format='png'):
        self.im.save(fh,format)

  def make_circle(self,draw, center, radius,colour=(0,255,0)):
    """ Cicle gradient is created by creating smaller and smaller cicles """
    (center_x, center_y) = center
    for i in range(0,radius):
      draw.ellipse(
        (center_x - radius + i,
         center_y - radius + i,
         center_x + radius - i,
         center_y + radius - i
        ),
        colour +(255 * i/(radius * 2),)
      )

  def add_circle(self, center, radius, colour):
    """ Adding a new cicle is a matter of creating a new one in a empty layer
    and merging it with the current one 

    XXX: Very heavy code, should actually only work on the data arrays, instead
    of doing all the magic with high-level images """

    im_new = Image.new("RGBA", self.im.size)
    draw = ImageDraw.Draw(im_new)
    self.make_circle(draw, center, radius, colour)
    
    data2 = np.array(im_new)
    
    # Add channels to make new images
    self.data = self.data + data2
    self.im = Image.fromarray(self.data)



class LatLonDeg():
  """ Helper class for coordinate conversions """
  def __init__(self,lat_deg, lon_deg):
    self.lat = lat_deg
    self.lon = lon_deg
  def __str__(self):
    return "%.5f,%.5f" % (self.lat, self.lon)

  def deg_per_pixel(self,other,pixel_max):
    return(LatLonDeg(abs(self.lat - other.lat) / pixel_max, abs(self.lon - other.lon) / pixel_max))



# Convertions of tile XYZ to WSG coordinates stolen from:
#    http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
# <stolen>
import math
def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return(xtile, ytile)

def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return(LatLonDeg(lat_deg,lon_deg))
# </stolen>


def boundbox_deg(x,y,z):
  """ Calculate the boundingbox for a image """
  return (num2deg(x,y,z), num2deg(x+1,y+1,z))



def make_tile(x,y,z,filter={},colour=(255,0,0)):
  """
  Crude attempt to generate tiles, by placing a gradient circle on a
  coordinate point.  Generate a larger tile  and make sure to plot related
  points first and then crop it to the required size (250x250).
  
  Many stuff NOT implemented yet, like:
  - Caching Images.
  - Conditional Filtering of Meting to allow display of sub-results.
  - Defining a extra level of transparency if you like to layer multiple tiles
    on top of each-other.
  - Color variation, allow the user to dynamically choose a the colour the
    points to be.
  - Advanced data plotting, like trying to guess the remainder points.
  """
  
  SIZE = 250
  
  nw_deg,se_deg = boundbox_deg(x,y,z)

  
  Picture = PyGamePicture
  resolution_deg = nw_deg.deg_per_pixel(se_deg, SIZE)
  # Converting LatLon to Meters is discussed here:
  #  http://stackoverflow.com/questions/3024404/transform-longitude-latitude-into-meters
  tile_height = float(40008000) / (2 ** z)
  meters_per_pixel = float(tile_height) / SIZE
  
  # Worst case scenario could a circle with 100% 'outside' our 250x250 range
  # also add data to the picture as circles are used
  border_pixels = 100 / meters_per_pixel / 2

  im = Picture("RGBA", (SIZE + border_pixels * 2,) * 2)

  nw_deg.lat += resolution_deg.lat * border_pixels
  nw_deg.lon -= resolution_deg.lon * border_pixels
  se_deg.lat -= resolution_deg.lat * border_pixels
  se_deg.lon += resolution_deg.lon * border_pixels
  
  lat_min = 999
  lon_min = 999
  lat_max = 0
  lon_max = 0

  filter.update({
    'latitude__lte' : nw_deg.lat,
    'latitude__gte' : se_deg.lat,
    'longitude__lte' : se_deg.lon,
    'longitude__gte' : nw_deg.lon
  })
  # TODO: This is currently hard-coded to display _all_ metingen
  metingen = Meting.objects.select_related().filter(**filter)

  # XXX: Signal is not normalized in the database making it unknown when a
  # signal is said to be 100% or when it is actually less, currently seems to
  # copy the raw reported values
  MAX_SIGNAL = 50
  # XXX: The radius relates to the zoom-level we are in, and should represent
  # a fixed distance, given the scale. Assume signal/distance to be lineair
  # such that signal 100% = 100m and 1% = 1m.
  # 
  # XXX: The relation is not lineair but from a more logeritmic scape, as we
  # are dealing with radio signals
  #
  MAX_RANGE = 100
  
  def dif(x,y):
    """ Return difference between two points """
    return max(x,y) - min(x,y)
  
  for meting in metingen:
    lat_min = min(lat_min, meting.latitude)
    lat_max = max(lat_max, meting.latitude)
    lon_min = min(lon_min, meting.longitude)
    lon_max = max(lon_max, meting.longitude)
    xcoord = int(dif(nw_deg.lon,meting.longitude) / (resolution_deg.lon))
    ycoord = int(dif(nw_deg.lat,meting.latitude) / (resolution_deg.lat))
    log.info(meting.accespoint.ssid, meting.latitude, meting.longitude, xcoord, ycoord)

    # TODO: Please note that this 'logic' technically does apply to WiFi signals,
    # if you are plotting from the 'source'. When plotting 'measurement' data you
    # get different patterns and properly need to start looking at techniques like:
    # Multilateration,Triangulation or Trilateration to recieve 'source' points.
    # 
    # Also you can treat all points as seperate and use techniques like
    # Multivariate interpolation to make the graphs. A nice overview at:
    #     http://en.wikipedia.org/wiki/Multivariate_interpolation
    #
    # One very intersting one to look at will be Inverse distance weighting
    # with examples like this: 
    # http://stackoverflow.com/questions/3104781/inverse-distance-weighted-idw-interpolation-with-python
    signal_normalized = MAX_RANGE - (MAX_SIGNAL - meting.signaal)
    im.add_circle((xcoord,ycoord),float(signal_normalized) / meters_per_pixel,colour, MAX_SIGNAL - meting.signaal)
  
  log.info("BoundingBox NW: %s" % nw_deg)
  log.info("BoundingBox SE: %s" % se_deg)
  log.info("")
  log.info("MetingBox   NW: %.5f,%.5f" % (lat_max, lon_min))
  log.info("MetingBox   SE: %.5f,%.5f" % (lat_min, lon_max))
  log.info("")
  log.info("Metingen Count: %s" % metingen.count())

  im.center_crop((SIZE,SIZE))
  return im


# Create your views here.
def serve_tile(request,zoom,x,y):
  filter = {}
  colour = (255,0,0)
  for key, value in request.GET.iteritems():
    if key == 'colour':
      colour = tuple(map(int,value.split(',')))
    else:
      filter[key] = value
  im = make_tile(int(x),int(y),int(zoom),filter=filter,colour=colour)
  response = HttpResponse(mimetype="image/png")
  im.write(response,'png')
  return response