generate-points.py

from random import randint
import numpy as np
import gdal
import osr
import psycopg2
import requests

def reverse_geocode_osm(coord):
    return requests.get(url = "http://nominatim.openstreetmap.org/reverse", params = dict(
        format='json',
        lon=coord[0],
        lat=coord[1],
        zoom='18',
        addressdetails='1'
    )).json()["address"]

conn = psycopg2.connect("dbname='hello' user='tlrobinson' host='localhost'")
cur = conn.cursor()

cur.execute("""CREATE TABLE IF NOT EXISTS people1 (
    longitude float,
    latitude float,
    country_code varchar,
    state varchar,
    city varchar,
    postcode varchar
);""")
conn.commit()

def insert(coord, address):
    params = (coord[1], coord[0], address.get("country_code"), address.get("state"), address.get("city"), address.get("postcode"))
    print params
    cur.execute("""INSERT INTO people1(longitude, latitude, country_code, state, city, postcode) VALUES (%s, %s, %s, %s, %s, %s)""", params)
    conn.commit()

ds = gdal.Open("gpw.tif")

width = ds.RasterXSize
height = ds.RasterYSize

source = osr.SpatialReference()
source.ImportFromWkt(ds.GetProjectionRef())

# create the new coordinate system
wgs84_wkt = """
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]"""
target = osr.SpatialReference()
target.ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(source, target)

#get the point to transform, pixel (0,0) in this case
gt = ds.GetGeoTransform()

def get_pixel_coord(index):
    return get_xy_pixel_coord(index % width, index / width)

def get_xy_pixel_coord(x, y):
    xx = gt[0] + x*gt[1] + y*gt[2]
    yy = gt[3] + x*gt[4] + y*gt[5]
    return transform.TransformPoint(xx, yy)

cached_filename = "population_cdf.npy"

try:
    cdf = np.load(cached_filename)
    print "loaded %s" % (cached_filename)
except Exception as e:
    print "creating %s" % (cached_filename)
    ds = gdal.Open("gpw.tif")
    raw = np.array(ds.GetRasterBand(1).ReadAsArray(), dtype=np.float64)
    cdf = np.cumsum(np.clip(raw, 0, None))
    print "saving %s" % (cached_filename)
    np.save(cached_filename, cdf)
    print "ok"

total = int(cdf[-1])

while True:
    value = randint(0, total)
    index = np.searchsorted(cdf, value)
    coord = get_pixel_coord(index)
    address = reverse_geocode_osm(coord)
    if address:
        insert(coord, address)