Generate dotmap tiles

dotmap.pde

/* I, Brandon Martin-Anderson, release this into the public domain or whatever. */


BufferedReader reader;
double ll, bb, rr, tt;

float A = 1000.0;

GlobalMercator proj = new GlobalMercator();

class PersonPoint {
  double x, y;
  String quadnode;

  PersonPoint(String row) {
    String[] fields = split(row, ",");
    this.x = Double.parseDouble(fields[0])/A;
    this.y = Double.parseDouble(fields[1])/A;
    this.quadnode = fields[2];
  }

  void draw(PGraphics pg) {
    pg.point((float)this.x, (float)this.y);
  }
}

ArrayList people;

float pointWeight(int level) {
  switch(level) {
  case 4: 
    return 0.05333;
  case 5: 
    return 0.08;
  case 6: 
    return 0.12;
  case 7: 
    return 0.18;
  case 8: 
    return 0.27;
  case 9: 
    return 0.405;
  case 10: 
    return 0.6075;
  case 11: 
    return 0.91125;
  case 12: 
    return 1.366875;
  case 13: 
    return 2.0503125;
  case 14: 
    return 3.07546875;
  case 15: 
    return 4.61320312;
  case 16:
    return 6.9198046;
  case 17:
    return 10.37970;
  default: 
    return 0.0;
  }
}

void setup() {
  size(512, 512, P2D);
  smooth();

  String[] zoomlevels = loadStrings("zoomlevel");

  for ( int i=0; i<zoomlevels.length; i++ ) {

    int level = int(zoomlevels[i]);

    println( "loading..." );
    reader = createReader("people.csv");
    try {
      String line;

      String quadkey = "";
      PGraphics pg = null;
      PVector tms_tile = null;
      
      int rown = 0;

      while (true) {
        line = reader.readLine();
        if (line==null || line.length()==0) {
          println( "file done" );
          break;
        }
        

        
        rown += 1;
        if( rown%100000==0 ){
          println( rown );
        }

        String[] fields = split(line, ",");
        float px = float(fields[0])/A;
        float py = float(fields[1])/A;
        String newQuadkey = fields[2].substring(0, level);
        
//          //only print out this quad:
//          if( !newQuadkey.substring(0,12).equals("032010110132") ){
//            continue;
//          }

        if ( !newQuadkey.equals( quadkey ) ) {
          //finish up the last tile
          if (pg!=null) {
            pg.endDraw();
            PVector gtile = proj.GoogleTile((int)tms_tile.x, (int)tms_tile.y, level);
            pg.save( "tiles/"+level+"/"+int(gtile.x)+"/"+int(gtile.y)+".png" );
            println( "done" );
          }

          quadkey = newQuadkey;
          
          PVector google_tile = proj.QuadKeyToTileXY( newQuadkey );
          tms_tile = proj.GoogleTile( (int)google_tile.x, (int)google_tile.y, level );

          println( level+" "+tms_tile.x+" "+tms_tile.y );

          pg = createGraphics(512, 512, P2D);
          pg.beginDraw();
          pg.smooth();

          PVector[] bounds = proj.TileBounds( (int)tms_tile.x, (int)tms_tile.y, level );

          float tile_ll = bounds[0].x/A;
          float tile_bb = bounds[0].y/A;
          float tile_rr = bounds[1].x/A;
          float tile_tt = bounds[1].y/A;

          double xscale = width/(tile_rr-tile_ll);
          double yscale = width/(tile_tt-tile_bb);
          float scale = min((float)xscale, (float)yscale);

          pg.scale(scale, -scale);
          pg.translate(-(float)tile_ll, -(float)tile_tt);

          pg.strokeWeight(pointWeight(level));

          pg.background(255);
        }

        pg.point(px, py);
      }

      if (pg!=null) {
        pg.endDraw();
        PVector gtile = proj.GoogleTile((int)tms_tile.x, (int)tms_tile.y, level);
        pg.save( "tiles/"+level+"/"+int(gtile.x)+"/"+int(gtile.y)+".png" );
        println( "done" );
      }
    } 
    catch (IOException e) {
      //e.printStackTrace();
    }
  }
}

void draw() {
}

gtile.pde

/* I, Brandon Martin-Anderson, release this into the public domain or whatever, unless I'm not allowed to, because it's actually a Java port of something I found here: http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/ */

class GlobalMercator {
  //  TMS Global Mercator Profile
  //	---------------------------
  //
  //	Functions necessary for generation of tiles in Spherical Mercator projection,
  //	EPSG:900913 (EPSG:gOOglE, Google Maps Global Mercator), EPSG:3785, OSGEO:41001.
  //
  //	Such tiles are compatible with Google Maps, Microsoft Virtual Earth, Yahoo Maps,
  //	UK Ordnance Survey OpenSpace API, ...
  //	and you can overlay them on top of base maps of those web mapping applications.
  //	
  //	Pixel and tile coordinates are in TMS notation (origin [0,0] in bottom-left).
  //
  //	What coordinate conversions do we need for TMS Global Mercator tiles::
  //
  //	     LatLon      <->       Meters      <->     Pixels    <->       Tile     
  //
  //	 WGS84 coordinates   Spherical Mercator  Pixels in pyramid  Tiles in pyramid
  //	     lat/lon            XY in metres     XY pixels Z zoom      XYZ from TMS 
  //	    EPSG:4326           EPSG:900913                                         
  //	     .----.              ---------               --                TMS      
  //	    /      \     <->     |       |     <->     /----/    <->      Google    
  //	    \      /             |       |           /--------/          QuadTree   
  //	     -----               ---------         /------------/                   
  //	   KML, public         WebMapService         Web Clients      TileMapService
  //
  //	What is the coordinate extent of Earth in EPSG:900913?
  //
  //	  [-20037508.342789244, -20037508.342789244, 20037508.342789244, 20037508.342789244]
  //	  Constant 20037508.342789244 comes from the circumference of the Earth in meters,
  //	  which is 40 thousand kilometers, the coordinate origin is in the middle of extent.
  //      In fact you can calculate the constant as: 2 * math.pi * 6378137 / 2.0
  //	  $ echo 180 85 | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:900913
  //	  Polar areas with abs(latitude) bigger then 85.05112878 are clipped off.
  //
  //	What are zoom level constants (pixels/meter) for pyramid with EPSG:900913?
  //
  //	  whole region is on top of pyramid (zoom=0) covered by 256x256 pixels tile,
  //	  every lower zoom level resolution is always divided by two
  //	  initialResolution = 20037508.342789244 * 2 / 256 = 156543.03392804062
  //
  //	What is the difference between TMS and Google Maps/QuadTree tile name convention?
  //
  //	  The tile raster itself is the same (equal extent, projection, pixel size),
  //	  there is just different identification of the same raster tile.
  //	  Tiles in TMS are counted from [0,0] in the bottom-left corner, id is XYZ.
  //	  Google placed the origin [0,0] to the top-left corner, reference is XYZ.
  //	  Microsoft is referencing tiles by a QuadTree name, defined on the website:
  //	  http://msdn2.microsoft.com/en-us/library/bb259689.aspx
  //
  //	The lat/lon coordinates are using WGS84 datum, yeh?
  //
  //	  Yes, all lat/lon we are mentioning should use WGS84 Geodetic Datum.
  //	  Well, the web clients like Google Maps are projecting those coordinates by
  //	  Spherical Mercator, so in fact lat/lon coordinates on sphere are treated as if
  //	  the were on the WGS84 ellipsoid.
  //	 
  //	  From MSDN documentation:
  //	  To simplify the calculations, we use the spherical form of projection, not
  //	  the ellipsoidal form. Since the projection is used only for map display,
  //	  and not for displaying numeric coordinates, we don't need the extra precision
  //	  of an ellipsoidal projection. The spherical projection causes approximately
  //	  0.33 percent scale distortion in the Y direction, which is not visually noticable.
  //
  //	How do I create a raster in EPSG:900913 and convert coordinates with PROJ.4?
  //
  //	  You can use standard GIS tools like gdalwarp, cs2cs or gdaltransform.
  //	  All of the tools supports -t_srs 'epsg:900913'.
  //
  //	  For other GIS programs check the exact definition of the projection:
  //	  More info at http://spatialreference.org/ref/user/google-projection/
  //	  The same projection is degined as EPSG:3785. WKT definition is in the official
  //	  EPSG database.
  //
  //	  Proj4 Text:
  //	    +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0
  //	    +k=1.0 +units=m +nadgrids=@null +no_defs
  //
  //	  Human readable WKT format of EPGS:900913:
  //	     PROJCS["Google Maps Global Mercator",
  //	         GEOGCS["WGS 84",
  //	             DATUM["WGS_1984",
  //	                 SPHEROID["WGS 84",6378137,298.2572235630016,
  //	                     AUTHORITY["EPSG","7030"]],
  //	                 AUTHORITY["EPSG","6326"]],
  //	             PRIMEM["Greenwich",0],
  //	             UNIT["degree",0.0174532925199433],
  //	             AUTHORITY["EPSG","4326"]],
  //	         PROJECTION["Mercator_1SP"],
  //	         PARAMETER["central_meridian",0],
  //	         PARAMETER["scale_factor",1],
  //	         PARAMETER["false_easting",0],
  //	         PARAMETER["false_northing",0],
  //	         UNIT["metre",1,
  //	             AUTHORITY["EPSG","9001"]]]

  int tileSize=256;
  float initialResolution;
  float originShift;

  GlobalMercator() {
    //Initialize the TMS Global Mercator pyramid"
    this.initialResolution = 2 * PI * 6378137 / this.tileSize;
    //156543.03392804062 for tileSize 256 pixels
    this.originShift = 2 * PI * 6378137 / 2.0;
    //20037508.342789244
  }

  PVector LatLonToMeters(float lat, float lon ) {
    //"Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913"

    float mx = lon * this.originShift / 180.0;
    float my = log( tan((90 + lat) * PI / 360.0 )) / (PI / 180.0);

    my = my * this.originShift / 180.0;
    return new PVector(mx, my);
  }

  PVector MetersToLatLon(float mx, float my ) {
    //"Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum"

    float lon = (mx / this.originShift) * 180.0;
    float lat = (my / this.originShift) * 180.0;

    lat = 180 / PI * (2 * atan( exp( lat * PI / 180.0)) - PI / 2.0);
    return new PVector(lat, lon);
  }

  PVector PixelsToMeters(float px, float py, int zoom) {
    //"Converts pixel coordinates in given zoom level of pyramid to EPSG:900913"

    float res = this.Resolution( zoom );
    float mx = px * res - this.originShift;
    float my = py * res - this.originShift;
    return new PVector(mx, my);
  }

  PVector MetersToPixels(float mx, float my, int zoom) {
    //"Converts EPSG:900913 to pyramid pixel coordinates in given zoom level"

    float res = this.Resolution( zoom );
    float px = (mx + this.originShift) / res;
    float py = (my + this.originShift) / res;
    return new PVector( px, py );
  }

  PVector PixelsToTile( float px, float py) {
    //"Returns a tile covering region in given pixel coordinates"

    int tx = int( ceil( px / float(this.tileSize) ) - 1 );
    int ty = int( ceil( py / float(this.tileSize) ) - 1 );
    return new PVector(tx, ty);
  }

  PVector PixelsToRaster(float px, float py, int zoom) {
    //"Move the origin of pixel coordinates to top-left corner"

    int mapSize = this.tileSize << zoom;
    return new PVector( px, mapSize - py );
  }

  PVector MetersToTile(float mx, float my, int zoom) {
    //"Returns tile for given mercator coordinates"

    PVector pp = this.MetersToPixels( mx, my, zoom);
    return this.PixelsToTile( pp.x, pp.y );
  }

  PVector[] TileBounds(int tx, int ty, int zoom) {
    //"Returns bounds of the given tile in EPSG:900913 coordinates"

    PVector ll = this.PixelsToMeters( tx*this.tileSize, ty*this.tileSize, zoom );
    PVector ur = this.PixelsToMeters( (tx+1)*this.tileSize, (ty+1)*this.tileSize, zoom );

    PVector[] ret = new PVector[2];
    ret[0]=ll;
    ret[1]=ur;
    return ret;
  }

  //	def TileLatLonBounds(self, tx, ty, zoom ):
  //		"Returns bounds of the given tile in latutude/longitude using WGS84 datum"
  //
  //		bounds = self.TileBounds( tx, ty, zoom)
  //		minLat, minLon = self.MetersToLatLon(bounds[0], bounds[1])
  //		maxLat, maxLon = self.MetersToLatLon(bounds[2], bounds[3])
  //		 
  //		return ( minLat, minLon, maxLat, maxLon )

  float Resolution( int zoom ) {
    //"Resolution (meters/pixel) for given zoom level (measured at Equator)"

    //# return (2 * math.pi * 6378137) / (self.tileSize * 2**zoom)
    return this.initialResolution / (pow(2, zoom));
  }

  //	def ZoomForPixelSize(self, pixelSize ):
  //		"Maximal scaledown zoom of the pyramid closest to the pixelSize."
  //		
  //		for i in range(30):
  //			if pixelSize > self.Resolution(i):
  //				return i-1 if i!=0 else 0 # We don't want to scale up
  //
  PVector GoogleTile(int tx, int ty, int zoom) {
    //"Converts TMS tile coordinates to Google Tile coordinates"

    // coordinate origin is moved from bottom-left to top-left corner of the extent
    return new PVector(tx, (pow(2, zoom) - 1) - ty);
  }
 
  //
  //	def QuadTree(self, tx, ty, zoom ):
  //		"Converts TMS tile coordinates to Microsoft QuadTree"
  //		
  //		quadKey = ""
  //		ty = (2**zoom - 1) - ty
  //		for i in range(zoom, 0, -1):
  //			digit = 0
  //			mask = 1 << (i-1)
  //			if (tx & mask) != 0:
  //				digit += 1
  //			if (ty & mask) != 0:
  //				digit += 2
  //			quadKey += str(digit)
  //			
  //		return quadKey
  
        PVector QuadKeyToTileXY(String quadKey)
        {
            int tileX;
            int tileY;
            int levelOfDetail;
          
            tileX = tileY = 0;
            levelOfDetail = quadKey.length();
            for (int i = levelOfDetail; i > 0; i--)
            {
                int mask = 1 << (i - 1);
                switch (quadKey.charAt(levelOfDetail - i))
                {
                    case '0':
                        break;

                    case '1':
                        tileX |= mask;
                        break;

                    case '2':
                        tileY |= mask;
                        break;

                    case '3':
                        tileX |= mask;
                        tileY |= mask;
                        break;

                }
            }
            
            return new PVector(tileX,tileY,levelOfDetail);
        }
}

Hi, me again! Can you post an example of your "zoomlevel" file?

Presumably his zoomlevel file looks something like:

4
5
6

loadStrings loads one line per file into a Strings array; level is used as int level = int(zoomlevels[i]);, and is used in 'substring' later to pull out the first n bytes of the quadkey.