1.表示するデータの取得
アメリカ統計局のサイトから人口データをダウンロードする。Shapefileも対応州をダウンロードする。基本データが公的機関で公開されているという点で、日本とは根本的に違う。
2.データの整形（コマンドラインによるshapefileの作成）
CSVCUTという非常にシンプルなツール（csvデータから必要部分を抜き出すだけの機能）と、ogr2ogrという地理データ（shapefile,TIGER,PostgreSQL, GML）の変換プログラムを使う。
dropdb census_demo
createdb census_demo -T template_postgis
# CSVCUTでデータをpopulation.csvに変換、PostgreSQLに取り込む
./csvcut -f 2,4,5 -o "," DEC_10_PL_GCTPL2.ST05/DEC_10_PL_GCTPL2.ST05.csv | tail +5 | cut -c 10-200 > population.csv
psql -d census_demo -c "create table population (id char(5), county varchar(100), population integer);"
psql -d census_demo -c "copy population from '${HOME}/src/making-maps-demo/population.csv' delimiters ',' CSV;"
# shapefileをPostgreSQLに取り込む
ogr2ogr -f PostgreSQL PG:dbname=census_demo tl_2010_17_county10/tl_2010_17_county10.shp -t_srs EPSG:900913 -nlt multipolygon -nln tl_2010_17_county10
# 合体して、「EPSG:4269」の単位系で書き出す
psql -d census_demo -f merge.sql
ogr2ogr -f "ESRI Shapefile" merged PG:dbname=census_demo -sql "select * from merged" -nln merged -t_srs EPSG:4269 -overwriteできたファイルをQGISで読み取ると以下のようになる。
3.地図タイルの自動生成
TileMillは、オープンソースの強力な地図作製ソフトウエア。Mapnikという別の強力なレンダリングライブラリ（OpenStreetMapに使われている）を利用している。
TileMill: Open Source Map Design from Development Seed on Vimeo.
TileMillでshapefileを読み込み、独自のCarto Languageという言語で地図の色分けをする。条件付き設定ができるので、データによって色を塗り分けたり、ズームレベルに応じてフォントサイズを調整する。
pythonのvirtualenvで仮想環境を作ってから、TileMillで作ったファイルをコピーして
cp path_to/TileMill/files/project/project_name/* tilemill
cartoというプロフラムで変換する
path_to/TileMill/bin/carto tilemill/project_name.mml > tilemill/project_name.xml
このxmlファイルに対して、OpenStreetMapのメンバーが作ったrender_tiles.pyという変換プログラムを修正したものを宛てる。レンダリングする範囲とズームレベル、使用するCPU数を指定する。
#!/usr/bin/python
from math import pi,cos,sin,log,exp,atan
from subprocess import call
import sys, os
from Queue import Queue
import mapnik2
import threading
import argparse
custom_fonts_dir = '/Library/Fonts/'
mapnik2.register_fonts(custom_fonts_dir)
DEG_TO_RAD = pi/180
RAD_TO_DEG = 180/pi
# Default number of rendering threads to spawn, should be roughly equal to number of CPU cores available
NUM_THREADS = 4
def minmax (a,b,c):
a = max(a,b)
a = min(a,c)
return a
class GoogleProjection:
def __init__(self,levels=18):
self.Bc = []
self.Cc = []
self.zc = []
self.Ac = []
c = 256
for d in range(0,levels):
e = c/2;
self.Bc.append(c/360.0)
self.Cc.append(c/(2 * pi))
self.zc.append((e,e))
self.Ac.append(c)
c *= 2
def fromLLtoPixel(self,ll,zoom):
d = self.zc[zoom]
e = round(d[0] + ll[0] * self.Bc[zoom])
f = minmax(sin(DEG_TO_RAD * ll[1]),-0.9999,0.9999)
g = round(d[1] + 0.5*log((1+f)/(1-f))*-self.Cc[zoom])
return (e,g)
def fromPixelToLL(self,px,zoom):
e = self.zc[zoom]
f = (px[0] - e[0])/self.Bc[zoom]
g = (px[1] - e[1])/-self.Cc[zoom]
h = RAD_TO_DEG * ( 2 * atan(exp(g)) - 0.5 * pi)
return (f,h)
class RenderThread:
def __init__(self, tile_dir, mapfile, q, printLock, maxZoom):
self.tile_dir = tile_dir
self.q = q
self.m = mapnik2.Map(256, 256)
self.printLock = printLock
# Load style XML
mapnik2.load_map(self.m, mapfile, True)
# Obtain Map projection
self.prj = mapnik2.Projection(self.m.srs)
# Projects between tile pixel co-ordinates and LatLong (EPSG:4326)
self.tileproj = GoogleProjection(maxZoom+1)
def render_tile(self, tile_uri, x, y, z):
# Calculate pixel positions of bottom-left & top-right
p0 = (x * 256, (y + 1) * 256)
p1 = ((x + 1) * 256, y * 256)
# Convert to LatLong (EPSG:4326)
l0 = self.tileproj.fromPixelToLL(p0, z);
l1 = self.tileproj.fromPixelToLL(p1, z);
# Convert to map projection (e.g. mercator co-ords EPSG:900913)
c0 = self.prj.forward(mapnik2.Coord(l0[0],l0[1]))
c1 = self.prj.forward(mapnik2.Coord(l1[0],l1[1]))
# Bounding box for the tile
if hasattr(mapnik2,'mapnik_version') and mapnik2.mapnik_version() >= 800:
bbox = mapnik2.Box2d(c0.x,c0.y, c1.x,c1.y)
else:
bbox = mapnik2.Envelope(c0.x,c0.y, c1.x,c1.y)
render_size = 256
self.m.resize(render_size, render_size)
self.m.zoom_to_box(bbox)
self.m.buffer_size = 128
# Render image with default Agg renderer
im = mapnik2.Image(render_size, render_size)
mapnik2.render(self.m, im)
im.save(tile_uri, 'png256')
def loop(self):
while True:
#Fetch a tile from the queue and render it
r = self.q.get()
if (r == None):
self.q.task_done()
break
else:
(name, tile_uri, x, y, z) = r
exists= ""
if os.path.isfile(tile_uri):
exists= "exists"
else:
self.render_tile(tile_uri, x, y, z)
bytes=os.stat(tile_uri)[6]
empty= ''
if bytes == 103:
empty = " Empty Tile "
self.printLock.acquire()
print name, ":", z, x, y, exists, empty
self.printLock.release()
self.q.task_done()
def render_tiles(bbox, mapfile, tile_dir, minZoom=1,maxZoom=18, name="unknown", num_threads=NUM_THREADS):
print "render_tiles(",bbox, mapfile, tile_dir, minZoom, maxZoom, name,")"
# Launch rendering threads
queue = Queue(32)
printLock = threading.Lock()
renderers = {}
for i in range(num_threads):
renderer = RenderThread(tile_dir, mapfile, queue, printLock, maxZoom)
render_thread = threading.Thread(target=renderer.loop)
render_thread.start()
#print "Started render thread %s" % render_thread.getName()
renderers[i] = render_thread
if not os.path.isdir(tile_dir):
os.mkdir(tile_dir)
gprj = GoogleProjection(maxZoom+1)
ll0 = (bbox[0],bbox[3])
ll1 = (bbox[2],bbox[1])
for z in range(minZoom,maxZoom + 1):
px0 = gprj.fromLLtoPixel(ll0,z)
px1 = gprj.fromLLtoPixel(ll1,z)
# check if we have directories in place
zoom = "%s" % z
if not os.path.isdir(tile_dir + zoom):
os.mkdir(tile_dir + zoom)
for x in range(int(px0[0]/256.0),int(px1[0]/256.0)+1):
# Validate x co-ordinate
if (x < 0) or (x >= 2**z):
continue
# check if we have directories in place
str_x = "%s" % x
if not os.path.isdir(tile_dir + zoom + '/' + str_x):
os.mkdir(tile_dir + zoom + '/' + str_x)
for y in range(int(px0[1]/256.0),int(px1[1]/256.0)+1):
# Validate x co-ordinate
if (y < 0) or (y >= 2**z):
continue
str_y = "%s" % y
tile_uri = tile_dir + zoom + '/' + str_x + '/' + str_y + '.png'
# Submit tile to be rendered into the queue
t = (name, tile_uri, x, y, z)
queue.put(t)
# Signal render threads to exit by sending empty request to queue
for i in range(num_threads):
queue.put(None)
# wait for pending rendering jobs to complete
queue.join()
for i in range(num_threads):
renderers[i].join()
if __name__ == "__main__":
#python render_tiles.py tilemill/wards.xml mayor-2011/.tiles/wards/ -89.03 41.07 -87.51 42.50 9 16 2
parser = argparse.ArgumentParser(description='Render your tiles.')
parser.add_argument('style_file', help="File containing Mapnik styles")
parser.add_argument('tile_dir', help="Destination directory for rendered tiles")
parser.add_argument('lat_1', type=float, help="Top-left corner latitude")
parser.add_argument('lon_1', type=float, help="Top-left corner longitude")
parser.add_argument('lat_2', type=float, help="Bottom-right corner latitude")
parser.add_argument('lon_2', type=float, help="Bottom-right corner longitude")
parser.add_argument('min_zoom', help="Minimum zoom level to render", type=int, default="9")
parser.add_argument('max_zoom', help="Maximum zoom level to render", type=int, default="12")
parser.add_argument('cores', help="Number of rendering threads to spawn, should be roughly equal to number of CPU cores available", type=int, default="2")
args = parser.parse_args()
bbox = (args.lon_1, args.lat_2, args.lon_2, args.lat_1)
render_tiles(bbox, args.style_file, args.tile_dir, args.min_zoom, args.max_zoom)
4.Google Mapに重ねる
プロは上手に書くものだ。
var map = null;
var center = new google.maps.LatLng(39.8, -89.2);
function fetch_tile(coord, zoom) {
return "./.tiles/" + zoom + "/" + coord.x + "/" + coord.y + ".png";
}
// --- begin geocoding stuff --
var geocoder = new google.maps.Geocoder();
var southwest_limit = new google.maps.LatLng(41.6, -88);
var northeast_limit = new google.maps.LatLng(42.05, -87.5);
var bounding_box = new google.maps.LatLngBounds(southwest_limit, northeast_limit);
var outside_il = false; // until prove true
var user_marker = null;
var has_zoomed = false;
var has_moved = false;
function geocode(query) {
if (typeof(query) == 'string') {
pattr = /\sil\s|\sillinois\s/gi;
match = query.match(pattr);
if (!match) {
query = query + ' IL';
}
gr = { 'address': query };
} else {
gr = { 'location': query };
}
geocoder.geocode(gr, handle_geocode);
}
function handle_geocode(results, status) {
alt_addresses(results);
lat = results[0].geometry.location.lat();
lng = results[0].geometry.location.lng();
normalized_address = results[0].formatted_address;
$('#query').val(normalized_address)
var zoom = (has_zoomed) ? map.zoom : 10;
process_location(lat, lng, zoom, true);
}
function alt_addresses(results) {
$('#alt-addresses').html('');
keep = new Array();
$.each(results, function(i,val) {
if (i==0) return; // skip the first result
for (var t in val.types) {
if (val.types[t] == 'street_address' || val.types[t] == 'intersection') {
keep.push(val.formatted_address);
break;
}
}
});
if (keep.length <= 1) {
$('#did-you-mean')
.addClass('disabled-link')
.unbind();
} else {
$('#did-you-mean')
.removeClass('disabled-link')
.click(function(e) {
e.stopPropagation();
toggle_alt_addresses();
});
$('#alt-addresses').append('<p>Did you mean...</p>');
for (var i in keep) {
$('#alt-addresses').append('<a href="javascript:geocode(\'' + keep[i] + '\');">' + keep[i] + '</a>');
}
}
}
function process_location(lat, lng, zoom, showMarker) {
var center = new google.maps.LatLng(lat, lng);
map.panTo(center);
if (zoom) {
map.setZoom(zoom);
}
if (showMarker instanceof Array) {
show_user_marker(showMarker[0],showMarker[1]); //?
} else if (showMarker) {
show_user_marker(lat, lng); //?
}
}
function show_user_marker(lat, lng) {
if (user_marker == null) {
user_marker = new google.maps.Marker();
user_marker.setMap(map);
}
user_marker.setPosition(new google.maps.LatLng(lat, lng));
}
function toggle_alt_addresses() {
alt_adds_div = $('#alt-addresses');
if (alt_adds_div.is(':hidden')) {
show_alt_addresses();
} else if (alt_adds_div.is(':visible')) {
hide_alt_addresses();
}
}
function show_alt_addresses() {
$('#alt-addresses').slideDown(250);
$('#did-you-mean').addClass('highlight');
}
function hide_alt_addresses() {
$('#alt-addresses').hide();
$('#did-you-mean.highlight').removeClass('highlight');
}
function parse_hash(s) {
try {
if (s[0] == "#") {
s = s.substr(1);
}
parts = s.split(",");
lat = parseFloat(parts[0]);
lng = parseFloat(parts[1]);
zoom = parseInt(parts[2]);
if (parts.length == 5) {
marker_location = [parts[3], parts[4]];
process_location(lat, lng, zoom, marker_location);
} else {
process_location(lat, lng, zoom, false);
}
} catch (e) { }
}
function make_hash() {
var parts = [map.center.lat(),map.center.lng(),map.zoom]
if (user_marker) {
parts.push(user_marker.position.lat());
parts.push(user_marker.position.lng());
}
return parts.join(",");
}
// --- end geocoding stuff ---
$(document).ready(function() {
census_demo_options = {
getTileUrl: fetch_tile,
tileSize: new google.maps.Size(256, 256),
isPng: true
}
census_demo = new google.maps.ImageMapType(census_demo_options);
map_options = {
minZoom: 7,
maxZoom: 10,
zoom: 7,
center: center,
mapTypeControl: false,
mapTypeId: "simple"
};
backdrop_styles = [
{
featureType: "administrative",
elementType: "labels",
stylers: [
{ lightness: 10 }
]
},{
featureType: "poi",
elementType: "labels",
stylers: [
{ visibility: "off" }
]
},{
featureType: "poi.park",
elementType: "geometry",
stylers: [
{ visibility: "off" }
]
},{
featureType: "road",
elementType: "geometry",
stylers: [
{ visibility: "simplified" },
{ saturation: -100 },
{ lightness: 0 }
]
},{
featureType: "road.arterial",
elementType: "labels",
stylers: [
{ gamma: 10 }
]
}
];
var initial_hash = window.location.hash;
simple = new google.maps.StyledMapType(backdrop_styles, { name: "Illinois population 2010" });
map = new google.maps.Map(document.getElementById("map_canvas"), map_options);
map.mapTypes.set("simple", simple);
map.overlayMapTypes.push(census_demo);
if (initial_hash) {
parse_hash(initial_hash);
}
$("#search").submit(function(){geocode($("#query").val());return false;});
$("#show-wards").click(function(){
if(this.checked){
map.overlayMapTypes.push(wards);
}
else {
map.overlayMapTypes.pop();
}
});
google.maps.event.addListener(map, 'zoom_changed', function() {
has_zoomed = true;
});
google.maps.event.addListener(map, 'bounds_changed', function() {
if (has_moved) {
window.location.hash = make_hash();
}
has_moved = true;
});
});