map generator see https://www.reddit.com/r/LosAngeles/comments/e767jc/an_artistic_map_of_la_freeways_using_30000000_gps/

process.py

#!/usr/bin/python

import sys
import png
import os
import numpy as np
import json
import math
import time
import random

from skimage.color import hsv2rgb
from itertools import chain

start = time.time()
max_for_line = int(os.environ.get('LMAX') or 50)
min_for_line = int(os.environ.get('LMIN') or 4)

pivot_ratio = math.sqrt(max_for_line - min_for_line)
pivot = pivot_ratio + min_for_line
ratio = float(os.environ.get('RATIO') or 0.3)

cx = 0
cy = 0

width = 9000 #* 3
height = 9000 #* 3

image = np.tile(np.tile(np.array([.5,.8,.0]), (width, 1)), (height, 1, 1))
sys.stderr.write("{:10} initialized\n".format( int(time.time() - start)) )

delta = [ 0.007, 0.01, 0.04 ]
dist = 0

delta_lng = .8 
delta_lat = delta_lng * (float(height)/width)

lng_low = -118.8
lng_high = lng_low + delta_lng

lat_low = 33.70
lat_high = lat_low + delta_lat

mult_lng = 360 / delta_lng * .999
mult_lat = 180 / delta_lat * .999
lng_comp = mult_lng * width / 360.0
lat_comp = mult_lat * height / 180.0

loc_last = False
mask = [x for x in delta]

def plotLineLow(x0,y0, x1,y1):
  dx = x1 - x0
  dy = y1 - y0
  yi = 1
  if dy < 0:
    yi = -1
    dy = -dy
  D = 2 * dy - dx
  y = y0

  frac = pivot/float(1 + dist - min_for_line)
  if frac > 1.5:
    mask[0] = frac * delta[0]
  else:
    mask[0] = -1/frac * delta[0]

  for x in range(x0, x1):
    if random.random() < (ratio * frac):
      image[y][x] = np.add(image[y][x], mask)
    if D > 0:
      y = y + yi
      D = D - 2 * dx
    D = D + 2 * dy

def plotLineHigh(x0,y0, x1,y1):
  dx = x1 - x0
  dy = y1 - y0
  xi = 1
  if dx < 0:
    xi = -1
    dx = -dx
  D = 2 * dx - dy
  x = x0

  frac = pivot/float(1 + dist - min_for_line)
  if frac > 1.5:
    mask[0] = frac * delta[0]
  else:
    mask[0] = -1/frac * delta[0]

  for y in range(y0, y1):
    if random.random() < (ratio * frac):
      image[y][x] = np.add(image[y][x], mask)
    if D > 0:
      x = x + xi
      D = D - 2*dy
    D = D + 2*dx

def convert(lat, lng):
  return int( lng_comp * (lng - lng_low)), int( lat_comp * (lat - lat_low))

with open(sys.argv[1], 'r') as json_file:
  points = json.load(json_file)

sys.stderr.write("{:10} json\n".format( int(time.time() - start) ) )

space = ''.join(['*'] * 10)
for booking in points[3000:]:
  loc_last = False
  for lat,lng in booking:

    cx += 1
    if cx % 100000 == 0:
      sys.stderr.write("{:17}".format("{}:{}:{} ".format(cx, cy * 1000 / cx, int(time.time() - start) )))
      if cx % 1000000 == 0:
        sys.stderr.write("\n")

    x, y = convert(lat, lng)
    
    if loc_last:
      x1, y1 = loc_last
      dy = abs(y1 - y)
      dx = abs(x1 - x)

      if dx + dy > min_for_line and dx + dy < max_for_line:
        x0, y0 = x, y

        # Tamil approx
        # 7/8dx + dy/2
        # We are scaling this back because we use it
        # for drawing
        #dist = dx * 7/8 + dy/2
        dist = math.sqrt( dx * dx + dy * dy)# , 0.5)

        #sys.stderr.write("{:10} {:17} {}\n".format(cx, dist, space[:int(dist/20)]))
        if dist > min_for_line:
          cy += 1
          if dy < dx:
            if x0 > x1:
              plotLineLow(x1, y1, x0, y0)
            else:
              plotLineLow(x0, y0, x1, y1)
            
          else:
            if y0 > y1:
              plotLineHigh(x1, y1, x0, y0)
            else:
              plotLineHigh(x0, y0, x1, y1)
    loc_last = [x, y]


sys.stderr.write("\n{:10} clipping\n".format( int(time.time() - start)) )
out = []
image = hsv2rgb(np.clip(image, 0, 1))

sys.stderr.write("{:10} scaling\n".format( int(time.time() - start)) )
for x in range(height):
  if (x+1) % 500 ==  0:
    sys.stderr.write("{:5} ".format(x+1))

  out.append(np.multiply(image[x], 255).astype(int).flatten().tolist())

sys.stderr.write("\n{:10} writing\n".format( int(time.time() - start) ) )

f = sys.stdout
w = png.Writer(width, height, greyscale=False)#, bitdepth=bitdepth)
w.write(f, reversed(out))
f.close()
sys.stderr.write("{:10} done\n".format( int(time.time() - start) ) )