| | 349 | |
| | 350 | ==== US Census ==== |
| | 351 | For the US, data is available down to block level: |
| | 352 | |
| | 353 | http://factfinder2.census.gov/faces/nav/jsf/pages/download_center.xhtml |
| | 354 | Dicennial Census |
| | 355 | 2010 SF1 100% Data |
| | 356 | |
| | 357 | Convert data from NAD83 to WGS84 (e.g. using qGIS) |
| | 358 | |
| | 359 | {{{ |
| | 360 | su postgres |
| | 361 | cd /data/Census2010BlockGroup |
| | 362 | shp2pgsql -s 4326 -I tl_2010_06037_bg10_WGS84.shp public.Census2010BlockGroup | psql -d gis |
| | 363 | psql |
| | 364 | \c gis |
| | 365 | ALTER TABLE census2010blockgroup ADD COLUMN population integer; |
| | 366 | ALTER TABLE census2010blockgroup ADD COLUMN population_density integer; |
| | 367 | \q |
| | 368 | exit |
| | 369 | |
| | 370 | w2p |
| | 371 | %autoindent |
| | 372 | pop_dict = {} |
| | 373 | input = os.path.join("/", "home", "data", "Census2010BlockGroup", "DEC_10_SF1_P1_with_ann.csv") |
| | 374 | inputFile = open(input, "r") |
| | 375 | header = 2 |
| | 376 | for line in inputFile: |
| | 377 | if header: |
| | 378 | header -= 1 |
| | 379 | continue |
| | 380 | parts = line.split(',', 6) |
| | 381 | geoid = parts[1] |
| | 382 | pop = int(parts[6].strip()) |
| | 383 | pop_dict[geoid] = pop |
| | 384 | |
| | 385 | inputFile.close() |
| | 386 | from __future__ import division |
| | 387 | db_string = "postgres://gis:GIS@localhost:5432/gis" |
| | 388 | db2 = DAL(db_string, migrate_enabled = False) |
| | 389 | table = db2.define_table("census2010blockgroup", Field("gid", "id"), Field("geoid10"), Field("aland10", "integer"), Field("population", "integer"), Field("population_density", "integer"), migrate=False) |
| | 390 | rows = db2().select(table.geoid10, table.aland10) |
| | 391 | for row in rows: |
| | 392 | data = {} |
| | 393 | geoid10 = row.geoid10 |
| | 394 | population = pop_dict.get(geoid10) |
| | 395 | if population: |
| | 396 | data["population"] = population |
| | 397 | aland10 = row.aland10 |
| | 398 | if aland10: |
| | 399 | area = aland10 / 2589988 |
| | 400 | population_density = population / area |
| | 401 | data["population_density"] = population_density |
| | 402 | db2(table.geoid10==geoid10).update(**data) |
| | 403 | |
| | 404 | db2.commit() |
| | 405 | }}} |
| | 406 | |
| | 407 | Can also have the Tract-level data for lower zoom (9-15) & create a Layer Group to serve the 2 layers with a style which shows only 1 level at each zoom |