Rails Against The Machine

This post has been superseded by (google heat maps in java)
I'm playing with the google maps API at the moment and thought it would be "fun" to try and do heat maps overlaid on google maps. This was hacked together in 2 days so its a bit rough round the edges.

Using Java for the server side stuff.

Part 1 javascript code for custom overlay.

follow the google api instructions to get a standard map working. Then create a new custom tile layer:

http://code.google.com/apis/maps/documentation/overlays.html#Tile_Overlays

function load() {
               if (GBrowserIsCompatible()) {
                   var myCopyright = new GCopyrightCollection("© ");
                   myCopyright.addCopyright(new GCopyright('Demo', new GLatLngBounds(new GLatLng(-90,-180), new GLatLng(90,180)),0,'©2007 Google'));
                  
                  
                   // Create the tile layer overlay and
                   // implement the three abstract methods                
                   var tilelayer = new GTileLayer(myCopyright);
                   tilelayer.getTileUrl = function(tile, zoom) {
                       return "http://localhost:8080/getTile?x="+tile.x+"&y="+tile.y+"&zoom="+zoom;
                   };
                   tilelayer.isPng = function() { return true;};
                   tilelayer.getOpacity = function() { return 0.75; }
                  
                   var myTileLayer = new GTileLayerOverlay(tilelayer);
                   var map = new GMap2(document.getElementById("map"));
                   map.addControl(new GSmallMapControl());
                   map.setCenter(new GLatLng(51.53608560178475, -0.06591796875), 14);//
                   map.addOverlay(myTileLayer);
               }
           }

Working out lat and long from google coordinates
Google uses the Mercator projection with units being pixels or tiles. This means that pixels map linearly to Longitude while Latitude is non linear equation given here

    private double pixel_to_lat(int pixel, int tile, int zoom) {
       double percent_down = (pixel + (height * tile))/(height * Math.pow(2, zoom));//work out what % of the way down we are
       double from_equator = 1f - 2f * percent_down;//but in standard mercator projection the equator is zero
       double y_mercator = from_equator * Math.PI;//and units run from + PI to -PI      
       double latitude = 2 * Math.atan(Math.exp(y_mercator)) - Math.PI / 2;
       return 180 * latitude / Math.PI;//convert to degrees
   }

   private double pixel_to_lng(int pixel, int tile, int zoom) {
       double circumference = width * Math.pow(2, zoom); //Number of tiles forming x axis
       double radius = circumference / (2 * Math.PI);
       double longitude = (pixel + width * tile) / radius;//in radians its just the ratio
       longitude = 180 * longitude / Math.PI;//convert to degrees
       return longitude - 180;//change center from date line to grenich
   }

so to check everything is working we draw squares with the lat and long written on them basically all you need to do is create a 255x255 png using the java buffered image thing.

    public static BufferedImage getTile(int xCoord, int yCoord, int zoom) {
       BufferedImage img = new BufferedImage(256, 256, BufferedImage.TYPE_INT_ARGB);
       Graphics2D g2 = img.createGraphics();
       g2.setStroke(new BasicStroke(1.5f));
       g2.setColor(Color.RED);
       g2.draw(new Rectangle2D.Double(0, 0, 256, 256));
       g2.setFont(new Font("Dialog", Font.BOLD, 25));
       g2.setColor(Color.RED);
       FontMetrics metrics = g2.getFontMetrics();
       DecimalFormat twoPlaces = new DecimalFormat("0.00");
       g2.drawString("x" + twoPlaces.format(tile_x_to_lng(xCoord, zoom)), 0, 0);
       g2.drawString("y" + twoPlaces.format(tile_y_to_lat(yCoord, zoom)), 0, metrics.getHeight());
       g2.drawString("z" + twoPlaces.format(zoom), 0, 2 * metrics.getHeight());
       return img;
   }

and to return this from our java server all we do is set the content type to png and write to the http response.

        int xCoord=Integer.parseInt(context.getRequest().getParameter("x"));
       int yCoord=Integer.parseInt(context.getRequest().getParameter("y"));
       int zoom=Integer.parseInt(context.getRequest().getParameter("zoom"));
      
       context.getResponse().setContentType("image/png");
       OutputStream output = context.getResponse().getOutputStream();
       ImageIO.write(ImageMaker.getTile(xCoord, yCoord, zoom), "png", output);
       //HeatMap h=new HeatMap();

So we can now confirm that it all worked:
add the following java script code http://code.google.com/apis/maps/documentation/examples/tile-detector.html

        GEvent.addListener(map,"click", function(overlay,latlng) {
         if (overlay) {
           // ignore if we click on the info window
           return;
         }
         var tileCoordinate = new GPoint();
         var tilePoint = new GPoint();
         var currentProjection = G_NORMAL_MAP.getProjection();
         tilePoint = currentProjection.fromLatLngToPixel(latlng, map.getZoom());
         tileCoordinate.x = Math.floor(tilePoint.x / 256);
         tileCoordinate.y = Math.floor(tilePoint.y / 256);
         var myHtml = "Latitude: " + latlng.lat() + "<br/>Longitude: " + latlng.lng() +
           "<br/>The Tile Coordinate is:<br/> x: " + tileCoordinate.x +
           "<br/> y: " + tileCoordinate.y + "<br/> at zoom level " + map.getZoom();
         map.openInfoWindow(latlng, myHtml);
       });

now hopefully the coordinates in the pop up should match the coordinates in the overlay. Something like:



Ok so now we know we can produce the correct map tiles on demand its time to do something more interesting like a heat map. Unfortunately the methods to draw individual pixels in java2d are really horrible.

 public static BufferedImage heatMap(int xCoord, int yCoord, int zoom, DataModel model) {
       int width = 256;
       int height = 256;
       byte[] pixels = generatePixels(width, height, model);
       //Put the results in a buffer so we can write it
       DataBuffer dbuf = new DataBufferByte(pixels, width * height, 0);
       //0x denotes hex so this 1111 i.e. take all the bits of the byte from 0 up to 15
       //and send them as the index to the indexed colour model
       int bitMasks[] = new int[]{(byte) 0xf};
       SampleModel sampleModel = new SinglePixelPackedSampleModel(DataBuffer.TYPE_BYTE, width, height, bitMasks);
       // Create a raster using the sample model and data buffer
       WritableRaster raster = Raster.createWritableRaster(sampleModel, dbuf, null);
       // Combine the color model and raster into a buffered image
       BufferedImage img = new BufferedImage(colorModel, raster, false, null);//new java.util.Hashtable());
       return img;
   }

   private static byte[] generatePixels(int w, int h, DataModel model) {
       byte[] pixels = new byte[w * h];
       int pIx = 0;
       for (int r = 0; r < h; r++) {
           for (int c = 0; c < w; c++) {
               pixels[pIx++] = (byte) (model.getValue(c, r));
           }
       }
       return pixels;
   }
   private static ColorModel colorModel = createColorModel();

   private static ColorModel createColorModel() {
// Generate 16-color model
       int n_col = 16;
       byte[] r = new byte[n_col];
       byte[] g = new byte[n_col];
       byte[] b = new byte[n_col];
       byte[] a = new byte[n_col]; //the alpha channel
       //Blue to red for first 8 levels
       for (int i = 0; i < 8; i++) {
           r[i] = (byte) (255 * i / 8);
           g[i] = (byte) 0;
           b[i] = (byte) (255 - 255 * i / 8);
           a[i] = (byte) (i * 255 / 16);
       }
       //red to yellow(255 255 0) for next 8 levels
       for (int i = 8; i < 16; i++) {
           r[i] = (byte) (255);
           g[i] = (byte) (255 * (i - 8) / 8);
           b[i] = (byte) 0;
           a[i] = (byte) (i * 255 / 16);
       }
       return new IndexColorModel(4, 16, r, g, b, a);
   }

Ok its horrible but the only really confusing bit is the bitmask. This is an array although in this case we only have one element. 0xf which means 1111 i.e use the first 4 bits to encode the colour index.

To use the class you just pass it a model which implements DataModel model where the datamodel interface is public double getValue(int x_pixels,int y_pixels); So if our model is something like:

    public double getValue(int x_pixel, int y_pixel) {
       //find the lat_long coordinates of the pixel
       double lon = pixel_to_lng(x_pixel,x_tile,z_tile);
       double lat = pixel_to_lat(y_pixel,y_tile,z_tile);
       //Now lets calculate the height from data lets just use a simple 2d gausian
       return gausian(lon, lat, -0.06591796875, 51.53608560178475);
   }

Then we get a pretty heat map of a symmetrical Gaussian.


Notice that the Gaussian seems longer in one direction than the other even though the distribution is symmetrical with respect to lat and long . This is because of the distorting effect of the projection and since our distance measure was just lat-lat_mean, lon-lon_mean rather than distance.

So lets try converting the distance measure to meters and drawing what should be a circle (with slight distortion due to the projection).

    private static double earthRadius = 6367500.0;

   public static double lat_lon_to_x_meters(double lat, double lon,double lon_origin) {
       return (Math.PI * earthRadius / 180.0) * (lon - lon_origin) * Math.cos(lat * (Math.PI / 180.0));
   }

   public static double lat_lon_to_y_meters(double lat, double lat_origin) {
       return (Math.PI * earthRadius / 180.0) * (lat - lat_origin);
   }

So lets draw a 100km radius circle and see what it looks like.


Next step kriging, Gaussian processes or SVMs to interpolate data or travel time maps like the ones Chris Lightfoot of http://www.mysociety.org/2007/more-travel-maps/ did. The big problem here is scalability we would like ideally to calculate them on the fly rather than as fixed one starting point images and since GPs require a matrix inversion there is obviously a bottle neck for large data sets.