ベクター地図タイル動的配信（Golang）

こんにちは。
下記のベクター地図タイル動的配信サーバを見つけたので動かしてみました。クライアント側の表示は MapboxGL を利用しました（そこの説明と同じです）[^1]。

[^1]: なお [Leaflet.VectorGrid](https://github.com/Leaflet/Leaflet.VectorGrid) 利用だと、現状はこの例の multi points データに対応できていないようです。

"[An example Go app for dynamically serving MapboxGL vector tiles](https://github.com/vicapow/go-vtile-example)" (GitHub)

* https://github.com/kkdd/go-vtile-example/tree/master （上記へ少し手を加えたもの）

```console
$ go run main.go
number of points = 91586
```

![go-vtile-example.jpg](https://qiita-image-store.s3.amazonaws.com/0/54617/d9a67881-13ce-08b5-a921-1475e8818304.jpeg)

### main.go

main.go に少し手を加えています。地理座標の内部表現（loc）はBraun投影した値です。

```math
\begin{align}
x &= \frac{\lambda}{2 \pi} \\
y &= \tan \frac{\phi}{2}
\end{align}
```

```main.go
package main

import (
	"os"
	"compress/gzip"
	"encoding/csv"
	"errors"
	"fmt"
	"io/ioutil"
	"log"
	"math"
	"net/http"
	"regexp"
	"strconv"
	"strings"
	"github.com/golang/protobuf/proto"
	"./vector_tile"
)

func cmdEnc(id uint32, count uint32) uint32 {
	return (id & 0x7) | (count << 3)
}

func moveTo(count uint32) uint32 {
	return cmdEnc(1, count)
}

func lineTo(count uint32) uint32 {
	return cmdEnc(2, count)
}

func closePath(count uint32) uint32 {
	return cmdEnc(7, count)
}

func paramEnc(value int32) int32 {
	return (value << 1) ^ (value >> 31)
}

func createTileWithPoints(points []XY, bounds XYZ) ([]byte, error) {
	layerName := "points"
	var pX, pY int32
	var layerVersion = vector_tile.Default_Tile_Layer_Version
	featureType := vector_tile.Tile_POINT
	var extent = vector_tile.Default_Tile_Layer_Extent
	fext := float64(extent);
	var geometry []uint32
	geometry = append(geometry, 0)  // npoints=0 (dummy)
	x, y := tileToBoundingBox(bounds)
	for _, point := range points {
		if point.x >= x[0] && point.x < x[1] && point.y >= y[0] && point.y < y[1] {
			p := locToTileXY(point, bounds)
			prevX := pX
			prevY := pY
			pX = int32(fext*p.x+0.5)
			pY = int32(fext*p.y+0.5)
			geometry = append(geometry, uint32(paramEnc(pX-prevX)))
			geometry = append(geometry, uint32(paramEnc(pY-prevY)))
		}
	}
	npoints := (uint32(len(geometry))-1)/2
	geometry[0] = moveTo(npoints)
	tile := &vector_tile.Tile{}
	tile.Layers = []*vector_tile.Tile_Layer{
		&vector_tile.Tile_Layer{
			Version: &layerVersion,
			Name:	&layerName,
			Extent:  &extent,
			Features: []*vector_tile.Tile_Feature{
				&vector_tile.Tile_Feature{
					Tags:	 []uint32{},
					Type:	 &featureType,
					Geometry: geometry,
				},
			},
		},
	}
	return proto.Marshal(tile)
}

// return loc: Braun projection
func lonLatToLoc(lonLat XY) (XY) {
	var loc XY
	loc.x = lonLat.x/360
	loc.y = math.Tan(lonLat.y/360 * math.Pi)  // Braun projection
	return loc
}

func locToLonLat(loc XY) (XY) {
	var lonLat XY
	lonLat.x = loc.x * 360
	lonLat.y = math.Atan(loc.y) * 360/math.Pi  // inverse Braun projection
	return lonLat
}

// relative position in a tile
func locToTileXY(loc XY, tile XYZ) (XY) {
	pos := loc
	pos.y = math.Log((1+pos.y)/(1-pos.y))/math.Pi/2  // web mercator
	pos.x = ( pos.x + 0.5) * tile.z - tile.x
	pos.y = (-pos.y + 0.5) * tile.z - tile.y
	return pos
}

func tileToLoc(tile XYZ) (XY) {
	var loc XY
	loc.x =   tile.x / tile.z - 0.5
	loc.y = -(tile.y / tile.z - 0.5)
	loc.y = 1 - 2/(math.Exp(loc.y*math.Pi*2)+1)  // inverse web mercator
	return loc
}

func tileToBoundingBox(tile XYZ) ([]float64, []float64) {
	upper := tileToLoc(tile)
	lower := tileToLoc(XYZ{x: tile.x, y: tile.y+1, z: tile.z})
	return []float64{upper.x, upper.x + 1/tile.z}, []float64{lower.y, upper.y}
}

const RE = 6378137.0  // GRS80
const FE = 1/298.257223563  // IS-GPS
const E2 = FE * (2 - FE)

//  geographic distance between two points
//  inputs: p = lonLatToLoc(lonLat1), q = lonLatToLoc(lonLat2)
func distance(p XY, q XY) (float64) {
	y2 := square((p.y + q.y) / 2)
	coslat := (1 - y2) / (1 + y2)
	w2 := 1 / (1 - E2 * (1 - coslat * coslat))
	dx := (p.x - q.x) * coslat
	dy := (p.y - q.y) * 2 / (1 + y2) * w2 * (1 - E2)
	return math.Sqrt(hypotSquared(dx, dy) * w2) * 2 * math.Pi * RE
}

func square(x float64) (float64) {
	return x * x
}

func hypotSquared(x float64, y float64) (float64) {
	return x * x + y * y
}

// Takes a string of the form `<z>/<x>/<y>` (for example, 1/2/3) and returns
// the individual uint32 values for x, y, and z if there was no error.
// Otherwise, err is set to a non `nil` value and x, y, z are set to 0.
func pathToTile(path string) (XYZ, error) {
	xyzReg := regexp.MustCompile("(?P<z>[0-9]+)/(?P<x>[0-9]+)/(?P<y>[0-9]+)")
	matches := xyzReg.FindStringSubmatch(path)
	if len(matches) == 0 {
		return XYZ{}, errors.New("Unable to parse path as tile")
	}
	x, err := strconv.ParseUint(matches[2], 10, 32)
	if err != nil {
		return XYZ{}, err
	}
	y, err := strconv.ParseUint(matches[3], 10, 32)
	if err != nil {
		return XYZ{}, err
	}
	z, err := strconv.ParseUint(matches[1], 10, 32)
	if err != nil {
		return XYZ{}, err
	}
	return XYZ{x: float64(x), y: float64(y), z: math.Pow(2, float64(z))}, nil
}

// A XYZ is a struct that holds tile's coordinates and zoom scale.
type XYZ struct {
	x float64
	y float64
	z float64
}

// A XY is a struct that holds a geographic location.
type XY struct {
	x float64
	y float64
}

// Tree a struct holder for tree information.
type Tree struct {
	lonlat XY
	species string
}

// trees.csv: TreeID,qLegalStatus,qSpecies,qAddress,SiteOrder,qSiteInfo,PlantType,qCaretaker,qCareAssistant,PlantDate,DBH,PlotSize,PermitNotes,XCoord,YCoord,Latitude,Longitude,Location
const SPECIES = 2
const LATITUDE = 15
const LONGITUDE = 16

func loadTrees() []Tree {
	content, err := ioutil.ReadFile("./trees.csv")
	if err != nil {
		log.Fatal(err)
	}
	r := csv.NewReader(strings.NewReader(string(content[:])))
	records, err := r.ReadAll()
	if err != nil {
		log.Fatal(err)
	}
	var trees []Tree
	for _, record := range records[1:] {
		species := record[SPECIES]
		lon, _ := strconv.ParseFloat(record[LONGITUDE], 64)
		lat, _ := strconv.ParseFloat(record[LATITUDE], 64)
		trees = append(trees, Tree{lonlat: XY{x: lon, y: lat}, species: species})
	}
	return trees
}

func main() {
	trees := loadTrees()
	points := make([]XY, len(trees), len(trees))
	for i, tree := range trees {
		points[i] = lonLatToLoc(tree.lonlat)
	}
	fmt.Println("number of points =", len(points))
	mux := http.NewServeMux()
	// Handle requests for urls of the form `/tiles/{z}/{x}/{y}` and returns
	// the vector tile for the even tile x, y, and z coordinates.
	tileBase := "/tiles/"
	mux.HandleFunc(tileBase, func(w http.ResponseWriter, r *http.Request) {
		log.Printf("url: %s", r.URL.Path)
		tile, err := pathToTile(r.URL.Path[len(tileBase):])
		if err != nil {
			http.Error(w, "Invalid tile url", 400)
			return
		}
		data, err := createTileWithPoints(points, tile)
		if err != nil {
			log.Fatal("error generating tile", err)
		}
		// All this APi to be requests from other domains.
		w.Header().Set("Content-Type", "application/x-protobuf")
		w.Header().Set("Access-Control-Allow-Origin", "*")
		w.Header().Set("Access-Control-Allow-Methods", "GET, POST, OPTIONS")
		w.Write(data)
	})
	log.Fatal(http.ListenAndServe(":8080", mux))
}
```