path: root/vodaviz.pde

//vodafone call visualisation
//
//Tim Redfern Jan 2012
//
//scale of wall 10.57x2.8m - some parts are obscured
//
//416 x 110 in. ~ 62400x16500 ~ 1029 mpx ~ 2945 MB ram to open
//
//CousinMarriageWorld.svg, nominally 940 × 470 pixels, file size: 1.99 MB
//World_map_(Miller_cylindrical_projection,_blank).svg‎, nominally 634 × 477 pixels, file size: 1.84 MB
//wikimedia
//
//http://www.vectortemplates.com/vector-world-map.php

//requirements - 
//must be able to transform points to a new projection
// maybe not through normal shape library
//must be able to find a random point within the shape
// 

//250112 so far so good, can iterate children.
//safest thing is to establish that I can place a point as being inside or outside of a shape
//(are they closed properly!)

//If we can find points in them then we can start the main task:
//connecting parts of the SVG with rows in the spreadsheet

//maybe next first I should look at projections (to keep them happy)
//I have a map with rectangular coords,
//is easy to enough to play with the projection

// public boolean contains(RGeomElem shp)
// public geomerative.RRectangle getBounds()
// --> this is kind of hidden. why? anyway, it should be possible to use it
//another option is to use the OUTLINE of the country (as a purely geometric thing).

/*
 http://en.wikipedia.org/wiki/3D_projection
 
 a- the 3D position of a point A that is to be projected.
 c- the 3D position of a point C representing the camera.
 O- The orientation of the camera (represented, for instance, by Tait–Bryan angles).
 e- the viewer's position relative to the display surface.[1]


http://williams.best.vwh.net/avform.htm#Crs
==good stuff

Intermediate points on a great circle
Here we find points (lat,lon) a given fraction of the distance (d) between them.
Suppose the starting point is (lat1,lon1) and the final point (lat2,lon2) and we want the point
a fraction f along the great circle route. f=0 is point 1. f=1 is point 2. The two points cannot
be antipodal ( i.e. lat1+lat2=0 and abs(lon1-lon2)=pi) because then the route is undefined.
The intermediate latitude and longitude is then given by:

        A=sin((1-f)*d)/sin(d)
        B=sin(f*d)/sin(d)
        x = A*cos(lat1)*cos(lon1) +  B*cos(lat2)*cos(lon2)
        y = A*cos(lat1)*sin(lon1) +  B*cos(lat2)*sin(lon2)
        z = A*sin(lat1)           +  B*sin(lat2)
        lat=atan2(z,sqrt(x^2+y^2))
        lon=atan2(y,x)
        
        //find 
        
        RPoint sp,ep,mp;
        float f=0.8; //fraction along path;
        float d=acos(sin(sp.y)*sin(ep.y)+cos(sp.y)*cos(ep.y)*cos(sp.x-ep.x));
        float A=sin((1-f)*d)/sin(d);
        float B=sin(f*d)/sin(d);
        float x = A*cos(lat1)*cos(lon1) +  B*cos(lat2)*cos(lon2);
        float y = A*cos(lat1)*sin(lon1) +  B*cos(lat2)*sin(lon2);
        float z = A*sin(lat1)           +  B*sin(lat2);
        mp.y=atan2(z,sqrt(x^2+y^2));
        mp.x=atan2(y,x);
        
        
        
        drawing - circles (fixed diameter), closed shape made of 2 bezier curves

 */

import processing.pdf.*;
import geomerative.*;

float log10 (float x) {
  return (log(x) / log(10));
}
class pointNormalise {
  //take pixel coords and turn into lat/lng radians
  float xo,xs,yo,ys;
  pointNormalise(float _xo,float _xs,float _yo,float _ys) {
    xo=_xo;
    xs=_xs;
    yo=_yo;
    ys=_ys;
  }
  RPoint alise(RPoint p) {
        //normalise
        float px=(((p.x-xo)/xs)-0.5)*PI*2; //26)/736);
        float py=(((p.y-yo)/ys)-0.5)*PI; //90)/390);
        
        return new RPoint(px,py);
  }
}

class pointTransform  {
  
  RPoint form(RPoint in) {
    //transform
    return new RPoint (((in.x*0.25)+(PI/2)),((in.y)*0.3)+(PI/2)); //front half of sphere
  }
  
}

class sphereMap {
        
    RPoint per(RPoint p,float _r) {
        

        //map to 3D sphere
        float r=getHeight()*_r;
        float x=r*cos(p.x)*(sin(p.y)) *2;
        float z=r*sin(p.x)*(sin(p.y))-985;
        float y=r*cos(p.y)+(getHeight()*0.5)-210;  
        
        //camera at 0,0,0 
        //screen plane at 0,0,100
        //(dx-ex)(ez/dz)
        //(dy-ey)(ez/dz)
        
        //rotate camera
        
        //transform into 2D plane @100

        return new RPoint(x*(100/z),y*(100/z));
    }
}

RPoint screenMapper(RPoint p) {
    p.x=((p.x*getWidth())/(PI*2))+(getWidth()/2);
    p.y=((p.y*getHeight())/PI)+(getHeight()/2);
    return p;
}

RPoint plerp(RPoint s,RPoint e,float a) {
  return new RPoint(lerp(s.x,e.x,a),lerp(s.y,e.y,a));
}

float GSphereDist(RPoint p1,RPoint p2) {
  return acos(sin(p1.y)*sin(p2.y)+cos(p1.y)*cos(p2.y)*cos(p1.x-p2.x));
}

RPoint GCircFract(RPoint sp,RPoint ep,float f) {
        float d=GSphereDist(sp,ep);
        float A=sin((1-f)*d)/sin(d);
        float B=sin(f*d)/sin(d);
        float x = A*cos(sp.y)*cos(sp.x) +  B*cos(ep.y)*cos(ep.x);
        float y = A*cos(sp.y)*sin(sp.x) +  B*cos(ep.y)*sin(ep.x);
        float z = A*sin(sp.y)           +  B*sin(ep.y);
        return new RPoint(atan2(y,x),atan2(z,sqrt(pow(x,2)+pow(y,2))));
}

RShape shp;
pointNormalise pnorm;
pointTransform ptrans;
sphereMap smap;

String mode;

csvloader data;
calldata calls;

bitmapcountry Ireland;

PImage lightmap;

void setup(){
  println("vodaviz v0.21");
  RG.init(this);

  //mode="PDF";
  if (mode=="PDF") size(832,220,PDF, "vodaviz_test_080212.pdf"); //P3D); //832,220); //nb pdf is 800x600
  else size(832,220); //,PDF, "testoutput.pdf"); //P3D); //832,220); //nb pdf is 800x600
  
  smooth();
  float m = millis();
  
  shp = RG.loadShape("countries_named_mercator.svg"); //test_drawing.svg"); //world_countries_outlines_split.svg");
  
  pnorm = new pointNormalise(18.279,746.302,129.314,374.293); 
  ptrans = new pointTransform();
  smap = new sphereMap();
  
  lightmap=loadImage("earthlights2_dmsp_big.jpg");
  
  
  RG.ignoreStyles();
  println("loaded svg in "+((millis()-m)*.001)+" seconds");
  Ireland=new bitmapcountry("Ireland",0,0,shp.children[0]);
  Ireland.analyse(18.279,746.302,129.314,374.293);
  data=new csvloader("calls.csv");
  calls=new calldata(data.data,shp);
 
  background(255,255,255);
  noFill();
  stroke(255);
  strokeWeight(.02);
  
  if (false) { //check worked example http://williams.best.vwh.net/avform.htm#Example
    RPoint LAX=new RPoint(2.066470,0.592539);
    RPoint JFK=new RPoint(1.287762,0.709186);
    println("LAX to JFK:"+GSphereDist(LAX,JFK)+" radians");
    RPoint LmJ=GCircFract(LAX,JFK,0.4);
    println("40% of LAX to JFK:"+LmJ.y+","+LmJ.x+" radians");
  }
  

  
}

int i=0;
int j=0;

void draw() {
 //if (i==0) image(lightmap,0,0,getWidth(),getHeight());
  
  //pick a random colour
  stroke(random(150)+10,random(150)+10,random(150)+10);

  if (true) //draw globe
  {
       if (true) { //draw countries
        beginShape();
        for (int k=0;k<calls.getcountry(i).outline.paths[0].commands.length;k++) {
          RPoint sp=calls.getcountry(i).outline.paths[0].commands[k].startPoint;
          RPoint dp=smap.per(ptrans.form(pnorm.alise(sp)),4);
          
          vertex(dp.x+(getWidth()/2),dp.y+(getHeight()/2)); //,z); //z);
        }
        endShape();
        //println("drawing "+calls.getcountry(i).name+": "+calls.getcountry(i).outline.paths[0].commands.length+" points");
       }
       
      
      if (true) {
        for (int j=0;j<calls.countries.get(i).calls*.001;j++) {
          RPoint s=calls.countries.get(i).getpoint();
          RPoint e=Ireland.getpoint();
           if (s.x>0&&e.x>0){ //point found
            RPoint sp=ptrans.form(pnorm.alise(s));
            RPoint ep=ptrans.form(pnorm.alise(e));
            RPoint Sp=smap.per(sp,4);
            RPoint Ep=smap.per(ep,4);
            RPoint Mp=smap.per(GCircFract(sp,ep,0.8),4.1);
              beginShape();
                line(Sp.x+(getWidth()/2),Sp.y+(getHeight()/2),Mp.x+(getWidth()/2),Mp.y+(getHeight()/2));
                line(Mp.x+(getWidth()/2),Mp.y+(getHeight()/2),Ep.x+(getWidth()/2),Ep.y+(getHeight()/2));
              endShape();
          }
        }
          //println("plotting "+calls.countries.get(i).name+": "+calls.countries.get(i).calls+" calls");
        }
        
  }
  else //draw 2D
  {
     if (true) { //draw countries
        beginShape();
        for (int k=0;k<calls.getcountry(i).outline.paths[0].commands.length;k++) {
          RPoint sp=calls.getcountry(i).outline.paths[0].commands[k].startPoint;
          RPoint dp=screenMapper(pnorm.alise(sp));
          
          vertex(dp.x,dp.y); //,z); //z);
        }
        endShape();
        //println("drawing "+calls.getcountry(i).name+": "+calls.getcountry(i).outline.paths[0].commands.length+" points");
      }
     if (true) { //draw lines
        for (int j=0;j<log10(calls.countries.get(i).calls)*100;j++) {
          RPoint s=calls.countries.get(i).getpoint();
          RPoint e=Ireland.getpoint();
           if (s.x>0&&e.x>0){ //point found
            RPoint Sp=screenMapper(pnorm.alise(s));
            RPoint Ep=screenMapper(pnorm.alise(e));
            RPoint Mp=plerp(Sp,Ep,0.75);
              beginShape();
                //line(Sp.x+(getWidth()/2),Sp.y+(getHeight()/2),Mp.x+(getWidth()/2),Mp.y+(getHeight()/2));
                //line(Mp.x+(getWidth()/2),Mp.y+(getHeight()/2),Ep.x+(getWidth()/2),Ep.y+(getHeight()/2));
                bezier(Sp.x, Sp.y, Sp.x, Sp.y-((getHeight()-Sp.y)*.1), Ep.x, Ep.y-((getHeight()-Ep.y)*.1), Ep.x, Ep.y);
              endShape();
          }
        }
          //println("plotting "+calls.countries.get(i).name+": "+calls.countries.get(i).calls+" calls");
        }
  }

       
  i++;

  //}
  if (i==calls.countries.size()-1) {
        println("finished");
       noLoop();
       if (mode=="PDF") exit();
  }
}

void mousePressed() {
  
}