pypixgrid.py

import os
import base64
import json
import sys
import argparse
import psycopg2
import psycopg2.extras
import mapbox_vector_tile
import sqlite3
from pprint import pprint
from math import sqrt
import zlib

def deflate(data, compresslevel=9):
    compress = zlib.compressobj(
            compresslevel,        # level: 0-9
            zlib.DEFLATED,        # method: must be DEFLATED
            16 + zlib.MAX_WBITS,  # window size in bits:
                                  #   -15..-8: negate, suppress header
                                  #   8..15: normal
                                  #   16..30: subtract 16, gzip header
            zlib.DEF_MEM_LEVEL,   # mem level: 1..8/9
            0                     # strategy:
                                  #   0 = Z_DEFAULT_STRATEGY
                                  #   1 = Z_FILTERED
                                  #   2 = Z_HUFFMAN_ONLY
                                  #   3 = Z_RLE
                                  #   4 = Z_FIXED
    )
    deflated = compress.compress(data)
    deflated += compress.flush()
    return deflated

def inflate(data):
    decompress = zlib.decompressobj(
            16 + zlib.MAX_WBITS  # see above
    )
    inflated = decompress.decompress(data)
    inflated += decompress.flush()
    return inflated

class PostGISProvider:
    def __init__(self, options):
        conn_string = "host='%s' dbname='%s' user='%s'" % (options['pg_connection']['host'], options['pg_connection']['dbname'], options['pg_connection']['user'])
        if options['pg_connection']['password']:
            conn_string += " password='%s'" % options['pg_connection']['password']
        conn = psycopg2.connect(conn_string)
	DEC2FLOAT = psycopg2.extensions.new_type(
	    psycopg2.extensions.DECIMAL.values,
	    'DEC2FLOAT',
	    lambda value, curs: float(value) if value is not None else None)
	psycopg2.extensions.register_type(DEC2FLOAT)
	self.conn = conn
        self.cursor = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
	self.cursor.execute(open("postgis_functions.sql", "r").read())
	self.conn.commit()

    def request(self, sql):
        # execute sql request over PostgreSQL connection
        self.cursor.execute(sql)
        return [dict(x) for x in self.cursor.fetchall()]

    def execute(self, sql):
        # execute sql request over PostgreSQL connection
        self.cursor.execute(sql)
	self.conn.commit()
	return 1



class MbTileWriter:
    def __init__(self,options,dbprovider):
	self.options = options["output"]
	self.config  = options 
	self.conn = sqlite3.connect(self.options["layername"]+'.mbtiles')
	self.cursor = self.conn.cursor()
	
	

	# creation de la table avec les metadonnees
	try:
		self.cursor.execute("CREATE TABLE metadata (name text, value text);")
	except:
		pprint("Un fichier mbtiles de meme nom existe deja.")
		return 0
	self.cursor.execute("INSERT INTO metadata VALUES ('name','{name}')".format(name=self.options["layername"]))
	self.cursor.execute("INSERT INTO metadata VALUES ('type','overlay')")
	if("version" in self.options):
		version = self.options["version"];
	else:
		version = 0
	self.cursor.execute("INSERT INTO metadata VALUES ('version','{version}')".format(version=version))
	if("description" in self.options):
		description = self.options["description"];
	else:
		description = "generated by pypixgrad"
	self.cursor.execute("INSERT INTO metadata VALUES ('description','{descr}')".format(descr=description))
	self.cursor.execute("INSERT INTO metadata VALUES ('format','{fo}')".format(fo=self.options["format"]))
	if("attribution" in self.options):
		attribution = self.options["attribution"];
		self.cursor.execute("INSERT INTO metadata VALUES ('description','{attr}')".format(attr=attribution))
	
	nbsop=len(options["scale_operations"])-1
	self.cursor.execute("INSERT INTO metadata VALUES ('maxzoom', {maz})".format(maz=options["scale_operations"][0][0]))
	self.cursor.execute("INSERT INTO metadata VALUES ('minzoom', {miz})".format(miz=options["scale_operations"][nbsop][0]))
	
	# noms de variables
	self.cursor.execute("INSERT INTO metadata VALUES ('temporal variables', {v})".format(v="'"+json.dumps(options["data_format"]["temporal_variables"])+"'"))
	if ("context_variables" in  options["data_format"]):
		self.cursor.execute("INSERT INTO metadata VALUES ('context variables', {v})".format(v="'"+json.dumps(options["data_format"]["context_variables"])+"'"))


	# recuperation de l'echantillonage temporel
	re = """select distinct(d.{time_column}) as time from {data_table} as d order by d.{time_column};
	""".format(time_column=config["data_format"]["time_column"],
		data_table=config["data_format"]["data_table"])
	timeindex=dbprovider.request(re)
	if(str(timeindex[0]["time"].__class__.__name__)=="datetime"):
		timeindex=map(lambda d: d["time"].strftime("%Y-%m-%dT%H:%M:%S"),timeindex)
	else:
		timeindex=map(lambda d: d["time"],timeindex)
	self.cursor.execute("INSERT INTO metadata VALUES ('time', {t})".format(t="'"+json.dumps(timeindex)+"'"))
	
	# recuperation des min/max de toutes les variables
	if ("temporal_variables" in  options["data_format"]):
		minmax=', '.join(["min("+v["name"]+") as min_"+v["name"]+", max("+v["name"]+") as max_"+v["name"]  for v in options["data_format"]["temporal_variables"]])
		rerange = "select " + minmax + " from " + config["data_format"]["data_table"]
		vtranges=dbprovider.request(rerange)
		self.cursor.execute("INSERT INTO metadata VALUES ('temporal variables ranges', {v})".format(v="'"+json.dumps(vtranges)+"'"))

	if ("context_variables" in  options["data_format"]):
		minmax=', '.join(["min("+v["name"]+") as min_"+v["name"]+", max("+v["name"]+") as max_"+v["name"]  for v in options["data_format"]["context_variables"]])
		rerange = "select " + minmax + " from " + config["data_format"]["geom_table"]
		vcranges=dbprovider.request(rerange) 
		self.cursor.execute("INSERT INTO metadata VALUES ('context variables ranges', {v})".format(v="'"+json.dumps(vcranges)+"'"))

	print("Metadata writed")
	# creation de la table des tuiles
	self.cursor.execute("CREATE TABLE tiles (zoom_level integer, tile_column integer, tile_row integer, tile_data blob);")
	self.conn.commit()

    def writerangeinmeta(self,rangeo,sc,vname):
	self.cursor.execute("INSERT INTO metadata VALUES ('{vname} variables ranges zoom {z}', {v})".format(v="'"+json.dumps(rangeo)+"'"),vname=vname,z=sc)


    def write(self, tile,x,y,z):
	sql = '''INSERT INTO tiles (zoom_level, tile_column, tile_row, tile_data) VALUES(?, ?, ?, ?);'''
	if(self.options["format"]=='pbf'):
		ct = MVTile(tile,self.options["layername"],self.config)
		yfliped = 2**z-1-y 
		self.cursor.execute(sql,[z,x,yfliped,sqlite3.Binary(deflate(mapbox_vector_tile.encode(ct.getContent())))])
    def commit(self):
	self.conn.commit()
		
class FileWriter:
    def __init__(self,options,dbprovider):
	self.config  = options
	self.options = options["output"]
	self.directory = self.options["directory"]
        if self.directory != '':
            if not os.path.exists(self.directory):
                os.makedirs(self.directory)


	# metadonnees
	nbsop=len(options["scale_operations"])-1
	metadata = {"maxzoom":options["scale_operations"][0][0],'minzoom':options["scale_operations"][nbsop][0]}


	re = """select distinct(d.{time_column}) as time from {data_table} as d order by d.{time_column};
	""".format(time_column=config["data_format"]["time_column"],
		data_table=config["data_format"]["data_table"])
	timeindex=dbprovider.request(re)
	if(str(timeindex[0]["time"].__class__.__name__)=="datetime"):
		timeindex=map(lambda d: d["time"].strftime("%Y-%m-%dT%H:%M:%S"),timeindex)
	else:
		timeindex=map(lambda d: d["time"],timeindex)
	metadata["time_range"]=timeindex;
	metadata["ranges"]={}
	self.metadata = metadata;
    def writerangeinmeta(self,rangeo,sc,vname):
	if sc in self.metadata["ranges"]:
		self.metadata["ranges"][sc][vname]=rangeo
	else:
		self.metadata["ranges"][sc]={vname:rangeo}
	

    def write(self, tile,x,y,z):
	self.z = str(z)
        self.zdir = str(z)
        if self.directory != '':
            self.zdir = self.directory + '/' + self.zdir
        if self.zdir != '':
            if not os.path.exists(self.zdir):
                os.makedirs(self.zdir)
        self.x = str(x)
        self.xdir = self.zdir + '/' + self.x
        if not os.path.exists(self.xdir):
            os.makedirs(self.xdir)
        self.y = str(y)
	
	if(self.options["format"]=='json'):
		ct = GeoJSONTile(tile)
		filename = self.xdir + '/' + self.y + '.json'
		with open(filename, 'w') as outfile:
	    		json.dump(ct.getContent(), outfile, sort_keys=True, indent=4, separators=(',', ': '))
	    		# pretty print
            		#json.dump(ct.getContent(), outfile, sort_keys=True, indent=4, separators=(',', ': '))
	if(self.options["format"]=='pbf'):
		ct = MVTile(tile,self.options["layername"],self.config)
		filename = self.xdir + '/' + self.y + '.pbf'
		with open(filename, 'wb') as outfile:
	    		# json.dump(self.content, outfile)
	    		# pretty print
            		outfile.write(mapbox_vector_tile.encode(ct.getContent()))
    def commit(self):
	filename = self.directory + "/metadata.json" 
	with open(filename, 'w') as outfile:	
		json.dump(self.metadata,outfile)
	print("Metadata writed")
	return 1
    
class GeoJSONTile:
    def __init__(self, data):
	self.content = {"type":"FeatureCollection"}
	features = []
	for o in data:
		lat = float(o["geometry"].split(' ')[0][6:])
		lng = float(o["geometry"].split(' ')[1][:-1])
		o.pop("geometry")
		o.pop("x")
		o.pop("y")
		o.pop("z")
		if(str(o["time"][0].__class__.__name__)=="datetime"):
			o["time"]=map(lambda d: d.strftime("%Y-%m-%dT%H:%M:%S"),o["time"])
		if(str(o["time"][0].__class__.__name__)=="time"):
			o["time"]=map(lambda d: d.strftime("%H:%M:%S"),o["time"])
		no = {"type" : "Feature", "geometry" : {"type" : "Point", "coordinates": [lat,lng]},"properties":o}
		features.append(no)
	self.content["features"]=features


    def getContent(self):
	return self.content;
        


class MVTile:
    def __init__(self, data, layername,config):
        self.content = {"name":layername}
	features = []
	for o in data:
		geom = o.pop("geometry")
		time = o.pop("time")
		po = []
		if(str(time[0].__class__.__name__)=="datetime"):
			time=map(lambda d: d.strftime("%Y-%m-%dT%H:%M:%S"),time)
		if ("temporal_variables" in  config["data_format"]):
			for v in config["data_format"]["temporal_variables"]:
				vals = o.pop(v["name"]) 
				for i in range(len(time)):
					po.append((v["name"]+"_"+str(time[i]),vals[i]))	
		if ("context_variables" in  config["data_format"]):
			for vc in config["data_format"]["context_variables"]:
				po.append((vc["name"],o.pop(vc["name"])))
		po.append(("area",o.pop("area")))
		po.append(("area_projected",o.pop("area_projected")))
		no = {"geometry" : geom,"properties":dict(po)}
		features.append(no)
	self.content["features"]=features

    def getContent(self):
	return self.content;
   


if __name__ == "__main__":

	try :
		with open(sys.argv[1]) as config_file:    
    			config = json.load(config_file)
	except :
		pprint("Veuillez fournir un fichier de configuration valide")
		sys.exit(0)
	provider = PostGISProvider(config)
	

	# test de la grille d'entree et extraction des params geometrique
	re = """select ST_AsText(g.{geom_column}) as geom, g.{row} as row, g.{col} as col, ST_SRID(g.{geom_column}) as srid from {geom_table} as g limit 2;
	""".format(geom_column=config["data_format"]["geom_column"],
			row = config["data_format"]["row_column"],
			col = config["data_format"]["col_column"],
			geom_table=config["data_format"]["geom_table"])
	grid_sample=provider.request(re)
	geo_ex=grid_sample[0]["geom"]
	coli = grid_sample[0]["col"]
	rowj = grid_sample[0]["row"]
	# la grille doit etre constituee de polygones	
	if(geo_ex[0:7]!='POLYGON'):
		pprint('Probleme de geometrie, verifier la configuration')
		sys.exit(0)
	coords=geo_ex[9:-2].split(",")
	# la grille doit etre constituee de rectangles	
	if(len(coords)!=5):
		pprint('Probleme de geometrie, les geometries doivent etre des carres. Verifier la configuration')
		sys.exit(0)
	p=map(lambda c: map(float,c.split(' ')),coords)
	config["data_format"]["grid_cell_size"]=p[2][0]-p[0][0]	

	xmin = min(map(lambda c: c[0],p))
	ymax = max(map(lambda c: c[1],p))
	if(p[0][0]!=xmin or p[0][1]!=ymax):
		pprint('Probleme de geometrie, le premier point de la geometrie doit etre le coin Nord-Ouest. Verifier la configuration')
		sys.exit(0)


	config["data_format"]["grid_origin"]=[p[0][0]-coli*config["data_format"]["grid_cell_size"],p[0][1]-rowj*config["data_format"]["grid_cell_size"]]

	config["data_format"]["grid_srid"]=grid_sample[0]["srid"]

	vsql      = ', '.join([v["aggregation"]+'('+v["name"]+') as '+v["name"] for v in config["data_format"]["temporal_variables"]])

	if ("context_variables" in  config["data_format"]):
		vcsql     = ','+', '.join([v["aggregation"]+'('+v["name"]+') as '+v["name"] for v in config["data_format"]["context_variables"]])
		vcnames   = ','+', '.join([v["name"] for v in config["data_format"]["context_variables"]])
	else:
		vcsql     = ''
		vcnames   = ''
	if(config["output"]["storage"]=="mbtiles"):
		writer   = MbTileWriter(config,provider)
	else:
		writer   = FileWriter(config,provider)

	if not ("nbquantiles" in config["output"]):
		config["output"]["nbquantiles"]=6

	# creation des tables aggregees	
	print("Aggregated grids creation")
	for i in range(len(config["scale_operations"])):
		so = config["scale_operations"][i]
		if(i==0):
			current_data_table=config["data_format"]["data_table"]
			current_geom_table=config["data_format"]["geom_table"]
			current_cell_size=config["data_format"]["grid_cell_size"]
		else :
			current_data_table="data_table_agg"+str(config["scale_operations"][i-1][0])
			current_geom_table="geom_table_agg"+str(config["scale_operations"][i-1][0])
			current_cell_size=current_cell_size*config["scale_operations"][i][1]
		

		data_table_sql = """create temp table data_table_agg{scale} as 
			select ceil(d.{row}/{agg_factor}) as {row}, floor(d.{col}/{agg_factor}) as {col}, d.{time} as {time}, {vsql} 
			from {data_table} as d group by ceil(d.{row}/{agg_factor}), floor(d.{col}/{agg_factor}), {time}; 
		""".format(
			scale=so[0],
			agg_factor = so[1],
			data_table = current_data_table,
			vsql = vsql,
			row = config["data_format"]["row_column"],
			col = config["data_format"]["col_column"],			
			time = config["data_format"]["time_column"])
		provider.execute(data_table_sql)
		
		geom_table_sql = """create temp table geom_table_agg{scale} as
			with newgrid as (select ceil(g.{row}/{agg_factor}) as {row}, floor(g.{col}/{agg_factor}) as {col} {vcsql} from {geom_table} as g 
			group by ceil(g.{row}/{agg_factor}), floor(g.{col}/{agg_factor}))
			select {row}, {col},
			ST_GeomFromText('Polygon(('||{xc0}||' '|| {yc0}||','||{xc0}||' '||{yc1}||','||{xc1}||' '||{yc1}||','||{xc1}||' '||{yc0}||','||{xc0}||' '||{yc0}||'))',{srid}) 
			as {geom_column} {vcnames} from newgrid; 
		""".format(
			vcsql=vcsql,
			vcnames=vcnames,
			scale=so[0],
			agg_factor = so[1],
			geom_table = current_geom_table,
			row = config["data_format"]["row_column"],
			col = config["data_format"]["col_column"],			
			geom_column = config["data_format"]["geom_column"],
			srid = config["data_format"]["grid_srid"],
			yc0  = config["data_format"]["row_column"]+'*'+str(current_cell_size)+'+'+str(config["data_format"]["grid_origin"][1]),
			yc1  = "("+config["data_format"]["row_column"]+'-1)*'+str(current_cell_size)+'+'+str(config["data_format"]["grid_origin"][1]),
			xc0  = config["data_format"]["col_column"]+'*'+str(current_cell_size)+'+'+str(config["data_format"]["grid_origin"][0]),
			xc1  = "("+config["data_format"]["col_column"]+'+1)*'+str(current_cell_size)+'+'+str(config["data_format"]["grid_origin"][0]),
			)
		provider.execute(geom_table_sql)
	print("Aggregated grids created")
	if("verbose" in config["output"] and config["output"]["verbose"]):

		for i in range(len(config["scale_operations"])):
			so = config["scale_operations"][i]
			pprint(["Statistiques des tuiles niveau :",so])
			stat_table_sql="""select ToTileX(g.{geom_column},{scale}) as X, ToTileY(g.{geom_column},{scale}) as Y, count(*), ST_Area(TileBBox({scale},cast(ToTileX(g.{geom_column},{scale}) as int),cast(ToTileY(g.{geom_column},{scale}) as int))) as area from geom_table_agg{scale} as g group by ToTileX(g.{geom_column},{scale}) , ToTileY(g.{geom_column},{scale});""".format(
					scale=so[0],
					geom_column = config["data_format"]["geom_column"]
					)
			stat_table = provider.request(stat_table_sql)
			pixpertiles = map(lambda r: r["count"], stat_table)
			pprint("carreaux par tuile, moyenne : {moy}, max : {max}".format(max=max(pixpertiles),moy=sum(pixpertiles) / float(len(pixpertiles))))
			tilesareas = map(lambda r: r["area"], stat_table)
			mt = sum(tilesareas) / float(len(tilesareas))
			area_sql="""select AVG(ST_Area(g.{geom_column})) as area from geom_table_agg{scale} as g ;""".format(
					scale=so[0],
					geom_column = config["data_format"]["geom_column"]
					)
			area_table=provider.request(area_sql)
			mc = area_table[0]["area"]
			rs = mc/mt
			spixel = int(sqrt(rs*256*256))
			pprint("surface moyenne des tuiles : {mt} (grid unit)^2, surface moyenne des carreaux : {mc} (grid unit)^2, ratio :{rs}".format(mt = mt,mc=mc,rs=rs))
			pprint("soit environ {s}x{s} pixels par carreaux".format(s = spixel))
			
	varraggsql = ', '.join(['array_agg(cast('+v["name"]+' as float) order by d.'+config["data_format"]["time_column"]+') as '+v["name"] for v in config["data_format"]["temporal_variables"]])
	if ("context_variables" in  config["data_format"]):
		vcnames   = ',' + ', '.join(['g.'+v["name"] for v in config["data_format"]["context_variables"]])
	else:
		vcnames   = ''

	# export des tuiles
	print("Tiles export")
	for i in range(len(config["scale_operations"])):
		so = config["scale_operations"][i]
		# requetes pour recuperer les donnees mises en forme pour l'export
		if(config["output"]["format"]=="json"):
			tiles_table_sql = """select ToTileX(g.{geom_column},{scale}) as X, ToTileY(g.{geom_column},{scale}) as Y, cast({scale} as int) as Z, 
				ST_AsText(ST_transform(ST_Centroid(g.{geom_column}),4326)) as geometry, ST_Area(g.{geom_column}) as area_projected, ST_Area((ST_Transform(g.{geom_column},4326))::geography) as area, array_agg(d.{time_column} order by d.{time_column}) as time, {varraggsql} {vcnames}
				from geom_table_agg{scale} as g, data_table_agg{scale} as d 
				where d.{row}=g.{row} and d.{col}=g.{col} 
				group by g.{col}, g.{row}, g.{geom_column} {vcnames} order by X, Y, Z;
			""".format(
				vcnames= vcnames,
				scale=so[0],
				row = config["data_format"]["row_column"],
				col = config["data_format"]["col_column"],			
				geom_column = config["data_format"]["geom_column"],
				time_column = config["data_format"]["time_column"],
				varraggsql = varraggsql
				)
		else:
			tiles_table_sql = """select ToTileX(g.{geom_column},{scale}) as X, ToTileY(g.{geom_column},{scale}) as Y, cast({scale} as int) as Z, 
				ST_AsText(mvtProject({geom_column},{scale})) as geometry, ST_Area(g.{geom_column}) as area_projected, ST_Area((ST_Transform(g.{geom_column},4326))::geography) as area, array_agg(d.{time_column} order by d.{time_column}) as time, {varraggsql} {vcnames}
				from geom_table_agg{scale} as g, data_table_agg{scale} as d 
				where d.{row}=g.{row} and d.{col}=g.{col} 
				group by g.{col}, g.{row}, g.{geom_column} {vcnames} order by X, Y, Z;
			""".format(
				vcnames= vcnames,
				scale=so[0],
				row = config["data_format"]["row_column"],
				col = config["data_format"]["col_column"],			
				geom_column = config["data_format"]["geom_column"],
				time_column = config["data_format"]["time_column"],
				varraggsql = varraggsql
				)
		
		
		tiles_table = provider.request(tiles_table_sql)
		xc = yc = zc = -1
		nbtiles = 0
		current_tile = []	
		for r in tiles_table:
			if (int(r["x"])!=xc or int(r["y"])!=yc or r["z"]!=zc):
	
				# on ecrit la tuile precedente
				if(nbtiles > 0):
					writer.write(current_tile,xc,yc,zc)
				nbtiles = nbtiles + 1;
				current_tile = []
			current_tile.append(r)
			xc=int(r["x"])
			yc=int(r["y"])
			zc=r["z"]	
		# ecriture de la derniere ligne
		writer.write(current_tile,xc,yc,zc)
		# recuperation des quantiles de toutes les variables
		for v in config["data_format"]["temporal_variables"]:
			quantiles="""WITH q AS (SELECT {v}, ntile({nbq}) over (order by {v}) AS quantile FROM data_table_agg{scale} where {v}>0)
			SELECT max({v}) as value, quantile as quantile FROM q GROUP BY quantile ORDER BY quantile""".format(scale=so[0],v=v["name"],nbq=config["output"]["nbquantiles"])
	
			vtranges=provider.request(quantiles)
			writer.writerangeinmeta(vtranges,so[0],v["name"])

		if ("context_variables" in  config["data_format"]):
			for v in config["data_format"]["context_variables"]:
				quantiles="""WITH q AS (SELECT {v}, ntile({nbq}) over (order by {v}) AS quantile FROM geom_table_agg{scale})
				SELECT max({v}) as value, quantile as quantile FROM q GROUP BY quantile ORDER BY quantile""".format(scale=so[0],v=v["name"],nbq=config["output"]["nbquantiles"])
	
				vcranges=provider.request(quantiles)
				writer.writerangeinmeta(vcranges,so[0],v["name"])

	print("Tiles exported")
	writer.commit()