generate svgs for any topo area in the world with python

rasterOutline.py

# download https://www.bodc.ac.uk/data/open_download/gebco/gebco_2020/geotiff/ and put into ./gebco/ folder
# fyi rasterio only supports up to python 2.7

import rasterio
import matplotlib.pyplot as plt
import math
import glob
import numpy as np
import datetime
import cv2
import time
import svgwrite

#fuji
# lat, lon = 35.360143833142075, 138.72906966596523

#florida
# lat, lon = -46.79134718528371, -71.50743733336654
lat, lon = 46.85221834809852, -121.76045166860439

#everest
lat, lon = 27.98392998449862, 86.92548767073704

lat,lon = 20.692537194215614, -156.38891942698237

for tile in glob.glob('gebco/*.tif'):
  dataset = rasterio.open(tile)

  py, px = dataset.index(lon, lat)
  if px > 0 and px < dataset.shape[0] and py > 0 and py < dataset.shape[1]  :
    band1 = dataset.read(1)

    width = 800
    height = 800

    mat = band1[py- height:py + height, px - width:px + width]

    print(mat)
    minval = np.min(mat)
    maxval = np.max(mat) 
    print(minval)
    print(maxval)
    
    valueRange = maxval - minval

    dwg = svgwrite.Drawing('test.svg', profile='full')

    # plt.matshow(mat)
    # plt.show()
    # bigmat = np.zeros()
    scale = 16
    scale = 2
    matUpsampled = cv2.resize(mat, (width *2 * scale, height*2*scale) ,interpolation = cv2.INTER_LANCZOS4)

    ticks = 30
    start = time.time()
    start = time.time()
    for tick in range(ticks):
        lineValue = tick * valueRange / ticks  + minval
        edgeImage = np.where(matUpsampled > lineValue, 0, 255)
    
        imgray = edgeImage.astype(np.int8)
        image, contours, hierarchy = cv2.findContours(imgray.astype(np.uint8),cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)

        for thisContour in contours:
            pathstring = "M " + str(thisContour[0][0][0]) + " " + str(thisContour[0][0][1]) + " "
            for i in range(1, len(thisContour)):
                coord = thisContour[i][0]
                pathstring += "L " + str(coord[0]) + " " + str(coord[1]) + " "

            pathstring += "Z"

            dwg.add(dwg.path( d = pathstring,
                stroke="#000",
                fill="none",
                stroke_width=1)
            )

            if tick% 2 == 0:
                dwg.add(dwg.path( d = pathstring ,
                    stroke="#000",
                    fill="none",
                    stroke_width=1))

            if tick% 4 == 0:
                dwg.add(dwg.path( d = pathstring ,
                    stroke="#000",
                    fill="none",
                    stroke_width=1))
                    
            if tick% 8 == 0:
                dwg.add(dwg.path( d = pathstring ,
                    stroke="#000",
                    fill="none",
                    stroke_width=1))
                

    lineHeight = 18
    duration = time.time() - start

    styleCSS ="font-size:16px;"
    line1 = 'GEBCO Compilation Group (2020) GEBCO 2020 Grid (doi:10.5285/a29c5465-b138-234d-e053-6c86abc040b9)'
    line2 ="lat, lon : {:.4f}, {:.4f}; extent: +/-{}x{}px; {} isolines extracted in {:.2f}s; range [{}m, {}m]".format(lat, lon, width, height, ticks, duration, minval, maxval)
    line3 = '{}'.format(datetime.datetime.now())

    dwg.add(dwg.text(line1, insert=(0, scale *height* 2 + lineHeight), style =styleCSS))
    dwg.add(dwg.text(line2, insert=(0, scale* height* 2 + 2 * lineHeight), style =styleCSS))
    dwg.add(dwg.text(line3, insert=(0, scale * height* 2 + 3 * lineHeight) , style =styleCSS))

    dwg.save()