---
# System prepended metadata

title: GIS - 地圖開發
tags: [前端]

---

---
tags: 前端
---
# GIS - 地圖開發
> GIS 地理資訊系統(Geographic Information System)
## 概念
向量(vector) | 網格(raster)
-|-
縮放不失真、儲存使用量小 | 渲染效率高，容易切割資源檔

![](http://www.seekacreative.co.nz/uploads/3/0/6/7/30679021/vector-raster_orig.jpg)

### 座標參考系統CRS:
* 台灣常用的：TWD97, WGS84, TWD67
* latitude(lat)：緯度
* longitude(lng、lon)：經度
* altitude：高度
* 投影法：麥卡托(Mercator)(方向準確)
* 大圓距離[Haversine](https://www.baeldung.com/cs/haversine-formula)
* 線上地圖服務WMS：openstreet、stamen、cartodb、esri、mapbox、google
* 檔案格式
    * shapefile：地理資訊系統始祖ESRI公司發明的
    * KML：谷歌推出的XML結構
    * geojson：開源通用JSON結構
        * Point | MultiPoint
        * LineSting | MultiLineString
        * Polygon | MultiPolygon
    ```json
    {
        type:'FeatureCollection',
        bbox:[lng1,lat1,lng2,lat2],
        features:[
            {
                type:'feature',
                geometry:{
                    type:'LineString',
                    coordinates:[
                        [lng,lat],
                        [lng,lat]
                    ]
                },
                properties: {}
            }
        ]
    }
    ```
## WMS資源
圖層|網址
-|-
特殊管制 | https://eghouse.hccg.gov.tw/arcgis/rest/services/Tiled3857/Nature_policy/MapServer/tile/{z}/{y}/{x}
水源水質 | https://eghouse.hccg.gov.tw/arcgis/rest/services/Tiled3857/Nature_water/MapServer/tile/{z}/{y}/{x}
地質敏感 | https://eghouse.hccg.gov.tw/arcgis/rest/services/Tiled3857/Nature_geo/MapServer/tile/{z}/{y}/{x}
環保與汙染 | https://eghouse.hccg.gov.tw/arcgis/rest/services/Tiled3857/Nature_environment/MapServer/tile/{z}/{y}/{x}
坡度示意圖 | https://eghouse.hccg.gov.tw/arcgis/rest/services/Tiled3857/SlopeTW3857_fix/MapServer/tile/{z}/{y}/{x}
國土利用 | https://wmts.nlsc.gov.tw/wmts/LUIMAP/default/EPSG:3857/{z}/{y}/{x}
公有土地 | https://wmts.nlsc.gov.tw/wmts/LAND_OPENDATA/default/EPSG:3857/{z}/{y}/{x}
學校 | https://wmts.nlsc.gov.tw/wmts/SCHOOL/default/EPSG:3857/{z}/{y}/{x}
村里界 | https://wmts.nlsc.gov.tw/wmts/Village/default/EPSG:3857/{z}/{y}/{x}
鄉市鎮界 | https://wmts.nlsc.gov.tw/wmts/TOWN/default/EPSG:3857/{z}/{y}/{x}
縣市界 | https://wmts.nlsc.gov.tw/wmts/CITY/default/EPSG:3857/{z}/{y}/{x}
OSM | https://{a-c}.tile.openstreetmap.org/{z}/{x}/{y}.png
通用 | https://wmts.nlsc.gov.tw/wmts/EMAP5/default/EPSG:3857/{z}/{y}/{x}
灰階 | https://wmts.nlsc.gov.tw/wmts/EMAP01/default/EPSG:3857/{z}/{y}/{x}
航照 | https://wmts.nlsc.gov.tw/wmts/PHOTO2/default/EPSG:3857/{z}/{y}/{x}
Google衛星 | https://mts1.google.com/vt/lyrs=s@186112443&hl=x-local&src=app&x={x}&y={y}&z={z}&s=Galile
Google地形 | https://mts1.google.com/vt/lyrs=p@186112443&hl=x-local&src=app&x={x}&y={y}&z={z}&s=Galile
ESRI衛星 | https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}
交通路網 | https://wmts.nlsc.gov.tw/wmts/EMAPX99/default/EPSG:3857/{z}/{y}/{x}


## Leaflet + 氣象地圖
* 可用python的folium產生網頁檔
* 開源近10年，擴充設計完善(熱力圖、叢集)，適用於視覺化
* API設計相當簡單，適合入門
* [Leaflet](https://leafletjs.com/)
* [Leaflet.markercluster](https://github.com/Leaflet/Leaflet.markercluster)
* [Leaflet.heat](https://github.com/Leaflet/Leaflet.heat)
* [topojson](https://github.com/topojson/topojson)

引入套件
```html
<link rel="stylesheet" href="https://unpkg.com/leaflet@1.7.1/dist/leaflet.css" />
<script src="https://unpkg.com/leaflet@1.7.1/dist/leaflet.js"></script>
<link rel="stylesheet" href="https://unpkg.com/leaflet.markercluster@1.4.1/dist/MarkerCluster.Default.css" />
<script src="https://unpkg.com/leaflet.markercluster@1.4.1/dist/leaflet.markercluster.js"></script>
<script src="https://www.unpkg.com/leaflet.heat@0.2.0/dist/leaflet-heat.js"></script>
<script src="https://unpkg.com/topojson@3.0.2/dist/topojson.min.js"></script>
```
配置畫面
```htmlembedded=
<body>
    <div id="map"></div>
</body>
<style>
    html, body, #map {
        margin: 0;
        width: 100%;
        height: 100%;
    }
</style>
```
初始化地圖
```javascript=
const map = L.map('map').setView([24, 121], 10);
L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', { attribution: '氣象地圖' }).addTo(map);
const markerCluster = L.markerClusterGroup().addTo(map);
```
撰寫Leaflet的[擴充套件](https://gist.github.com/rclark/5779673)
```javascript=
L.TopoJSON = L.GeoJSON.extend({
    addData(jsonData) {
        if (jsonData.type === "Topology") {
            for (let key in jsonData.objects) {
                let geojson = topojson.feature(jsonData, jsonData.objects[key]);
                L.GeoJSON.prototype.addData.call(this, geojson);
            }
        }
        else L.GeoJSON.prototype.addData.call(this, jsonData);
    }
});
L.topoJson = (data, options)=>new L.TopoJSON(data, options);
```
查詢氣象API的[中英文名對照表](https://opendata.cwb.gov.tw/opendatadoc/DIV2/A0001-001.pdf)
選取多行`ctrl+alt+↓`、選取相同字`ctrl+d`
```javascript=
const ZH_TW = {
    lat: '緯度',
    lon: '經度',
    locationName: '測站名稱',
    stationId: '測站ID',
    obsTime: '觀測時間',
    CITY: '縣市',
    TOWN: '地區',
    ELE: '高度',
    WDIR: '風向',
    WDSD: '風速',
    TEMP: '溫度',
    HUMD: '相對濕度',
    PRES: '測站氣壓',
    H_24R: '日累積雨量',
    H_FX: '小時最大陣風風速',
    H_XD: '小時最大陣風風向',
    H_FXT: '小時最大陣風時間',
    D_TX: '本日最高溫',
    D_TXT: '本日最高溫發生時間',
    D_TN: '本日最低溫',
    D_TNT: '本日最低溫發生時間',
    CITY: '縣市',
    CITY_SN: '縣市編號',
    TOWN: '鄉鎮',
    TOWN_SN: '鄉鎮編號',
}
const NAME = ['WDIR', 'WDSD', 'TEMP', 'HUMD', 'PRES', 'H_24R', 'H_FX', 'H_XD', 'D_TX', 'D_TN']
```
取得氣象API
```javascript=
async function getWeather(markerCluster) {
    let json = await fetch('https://opendata.cwb.gov.tw/api/v1/rest/datastore/O-A0001-001?Authorization=rdec-key-123-45678-011121314').then(res => res.json())
    let { success, result: { fields }, records: { location } } = json
    return location.map(loc => {
        let { lat, lon, locationName, parameter, stationId, time: { obsTime }, weatherElement } = loc
        let addr = Object.fromEntries(parameter.map(x => [x.parameterName, x.parameterValue]))
        let weather = Object.fromEntries(weatherElement.map(x => [x.elementName, x.elementValue]))
        let data = { lat, lon, locationName, stationId, ...addr, ...weather, obsTime }
        let table = document.createElement('table')
        table.className = 'table'
        for (let key in data) {
            let row = table.insertRow()
            row.insertCell().textContent = ZH_TW[key]
            row.insertCell().textContent = data[key] == -99 ? '無資料' : data[key]
        }
        L.marker([lat, lon]).bindPopup(table).addTo(markerCluster)
        return data
    })
}
```
分層設色
```javascript=
function getChoroplenth(topoJson = {}, colorMap = {}, popupMap = {}) {
    let config = {
        style: (feature) => ({
            color: "#000",
            opacity: 1,
            weight: 1,
            fillColor: colorMap[feature.properties.name],
            fillOpacity: 0.8
        })
        , onEachFeature: (feature, layer) => layer.bindPopup(feature.properties.name + ':' + popupMap[feature.properties.name])
    }
    return L.topoJson(topoJson, config).addTo(map)
}
```
工具函式
```javascript=
// 色票轉換 https://stackoverflow.com/questions/36721830/convert-hsl-to-rgb-and-hex
function hslToHex(h, s, l) {
    l /= 100;
    const a = s * Math.min(l, 1 - l) / 100;
    const f = n => {
        const k = (n + h / 30) % 12;
        const color = l - a * Math.max(Math.min(k - 3, 9 - k, 1), -1);
        return Math.round(255 * color).toString(16).padStart(2, '0');
    };
    return `#${f(0)}${f(8)}${f(4)}`;
}
function groupBy(rows = [], key = '') {
    let group = {}
    for (let row of rows) {
        let rKey = row[key]
        if (!(rKey in group)) group[rKey] = []
        group[rKey].push(row)
    }
    return group
}
function mean(arr = []) {
    return arr.reduce((acc, cur) => acc + cur, 0) / arr.length
}
function normalize(arr = []) {
    let min = Math.min(...arr)
    let max = Math.max(...arr)
    return arr.map(x => (x - min) / (max - min))
}
```
主程式
```javascript=
async function mounted() {
    let topo = await fetch('https://raw.githubusercontent.com/g0v/twgeojson/master/json/twCounty2010merge.topo.json').then(res => res.text())
    topo = JSON.parse(topo.replace(/台/g, '臺').replace(/桃園縣/g, '桃園市'))
    let rows = await getWeather(markerCluster);
    const heat = L.heatLayer(rows.map(row => [row.lat, row.lon])).addTo(map);
    // 定義控制項
    const checkboxes = { '標記叢集': markerCluster, '熱力圖': heat }
    const radios = { '無分層設色': L.layerGroup() }

    // 群組化
    let cityGroup = groupBy(rows, 'CITY')
    for (let name of NAME) {
        let meanMap = {}
        let colorMap = {}
        let popupMap = {}
        for (let key in cityGroup) {
            let arr = cityGroup[key].map(row => parseFloat(row[name])).filter(x => x != -99 && !isNaN(x))
            meanMap[key] = mean(arr)
        }
        let means = Object.values(meanMap)
        let norms = normalize(means)
        Object.keys(meanMap).map((key, i) => {
            colorMap[key] = hslToHex(norms[i] * 360, 50, 50)
            popupMap[key] = means[i].toFixed(2)
        })
        radios[ZH_TW[name]] = getChoroplenth(topo, colorMap, popupMap);
    }
    map.fitBounds(Object.values(radios)[1].getBounds());
    L.control.layers(radios, checkboxes).addTo(map);
}
mounted()
```

## Google Map
> 谷歌地圖開發者服務，適合用於開發多功能資訊整合
* 地圖JS
* 街景
* 商家景點
* 地理編碼(geocoding)
* 路徑規劃(direction)

引入套件
```htmlembedded=
<div id="map"></div>
<script src="https://maps.googleapis.com/maps/api/js?key={API_KEY}&callback=initMap"></script>
<script src="https://developers.google.com/maps/documentation/javascript/examples/markerclusterer/markerclusterer.js"></script>
```
基本圖層
```javascript=
var map = new google.maps.Map(document.getElementById('map'),{center:{lat,lng},zoom:8,mapTypeId: 'terrain'})
// 標記叢集
var marker = new google.maps.Marker({position: {lat,lng}, map}) 
var imagePath = 'https://developers.google.com/maps/documentation/javascript/examples/markerclusterer/m'
var markerCluster = new MarkerClusterer(map, [marker], {imagePath});
// 熱力圖
new google.maps.visualization.HeatmapLayer({data: [{lat,lng}],dissipating:false,map: map})
// 匯入資料
map.data.addGeoJson(GeoJson)
// 事件監聽
map.addEventListener('click',e=>{e.latlng.lat()})

new google.maps.Polyline({
    path: flightPlanCoordinates,
    geodesic: true,
    strokeColor: '#FF0000',
    strokeOpacity: 1.0,
    strokeWeight: 2
  }).setMap(map) //加上setMap(null)移除、加上getLength取得長度(Polygon)

new google.maps.Rectangle({
  bounds: {north: 33.685, south: 33.671, east: -116.234,west: -116.251},
  editable: true
}).setMap(map);
```
路徑規劃
```javascript=
var request = {
    origin: {lat,lng}, //傳入關鍵字(用+號分隔)、經緯度、place_id
    destination: end,
    travelMode: 'DRIVING' //DRIVING, WALKING, BICYCLING, TRANSIT
}
var directionsService = new google.maps.DirectionsService();
var directionsDisplay = new google.maps.DirectionsRenderer();
directionsDisplay.setMap(map)
directionsService.route(request, function(result, status) {
    if (status == 'OK') { // 另一個狀態為"ZERO_RESULTS"
      directionsDisplay.setDirections(result);
      result[i].legs[j] // {start_address, end_address, distance{text}}
      var distance = result[0].legs.reduce((acc,cur)=>acc+cur.distance.value,0)
    }
})
```

## OpenLaysers
* 收錄了完整地圖功能模組，適合開發地圖編輯軟體
* 資料結構高達8層，開發相對不易

```graphviz
digraph{
    map->layers;
    map->controls;
    map->interactions;
    map->overlays;
    map->view;
    view->projection, center, zoom, rotation
    overlays->element, offset, positioning
    layers->source;
    source->raster;
    source->cluster;
    source->vector;
    cluster->vector;
    raster->XYZ;
    vector->feature;
    feature->style;
    feature->geom;
    geom->coord;
    geom->type;
    type->Point, LineString, MultiLineString, Polygon, MultiPolygon
}
```

## 手刻地圖 - 核心數學
```javascript=
const tilePixel = { w: 256, h: 256 }
const tileSize = 6378137 * Math.PI
const World = {
    bbox: [-tileSize, -tileSize, tileSize, tileSize],
    w: 2 * tileSize,
    h: 2 * tileSize,
}
const View = {
    bbox: [],
    center: { x: 13437548.485305637, y: 2772337.9239074644 },
    w: this.canvas.clientWidth,
    h: this.canvas.clientHeight,
    zoom: 9,
    minZoom: 0,
    maxZoom: 20,
}
class Point{
    toLngLat(){
        var d = 180 / Math.PI
        return [
            point[0] * d / 6378137,
            (2 * Math.atan(Math.exp(point[1] / 6378137)) - (Math.PI / 2)) * d
        ]
    }
    getArea(aoa) {
        var radius = 6371008.8, area = 0
        let x1 = aoa[aoa.length - 1][0]
        let y1 = aoa[aoa.length - 1][1]
        for (let i = 0; i < aoa.length; i++) {
            let x2 = aoa[i][0]
            let y2 = aoa[i][1]
            area += (x2 - x1) * Math.PI / 180 * (2 + Math.sin(y1 * Math.PI / 180) + Math.sin(y2 * Math.PI / 180))
            x1 = x2
            y1 = y2
        }
        return Math.abs(area * radius * radius / 2.0)
    }
    client2point(client) {
        var view = this.view,
            W = view.bbox[2] - view.bbox[0],
            H = view.bbox[1] - view.bbox[3]
        return [
            view.bbox[0] + client[0] * W / view.w,
            view.bbox[3] + client[1] * H / view.h
        ]
    }
    point2client(point) {
        var view = this.view,
            W = view.bbox[2] - view.bbox[0],
            H = view.bbox[1] - view.bbox[3]
        return [
            (point[0] - view.bbox[0]) / W * view.w,
            (point[1] - view.bbox[3]) / H * view.h
        ]
    }
}
class LngLat{
    toPoint(lnglat) {
        // SphericalMercator
        var d = Math.PI / 180,
            max = 85.0511287798,
            lat = Math.max(Math.min(max, lnglat[1]), -max),
            sin = Math.sin(lat * d);
        return [
            6378137 * lnglat[0] * d,
            6378137 * Math.log((1 + sin) / (1 - sin)) / 2
        ]
    }
    getLength(c1, c2){
        var r = 6371008.8,
            lat1 = c1[1] * Math.PI / 180,
            lat2 = c2[1] * Math.PI / 180,
            dlat = (lat2 - lat1) / 2,
            dlon = (c2[0] - c1[0]) * Math.PI / 180 / 2,
            a = Math.sin(dlat) * Math.sin(dlat) +
                Math.cos(lat1) * Math.cos(lat2) *
                Math.sin(dlon) * Math.sin(dlon)
        return 2 * r * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a))
    }
}
class Bounds{
    overlaps(bounds){
        let x = (bounds.max.x > this.min.x) && (bounds.min.x < this.max.x)
        let y = (bounds.max.y > this.min.y) && (bounds.min.y < this.max.y)
        return x&&y
    }
    contains(point){
        let x = (point.x > this.min.x) && (point.x < this.max.x)
        let y = (point.y > this.min.y) && (point.y < this.max.y)
        return x&&y
    }
}
function getBbox(aoa) {
    return [
        Math.min(...aoa.map(c => c[0])),
        Math.min(...aoa.map(c => c[1])),
        Math.max(...aoa.map(c => c[0])),
        Math.max(...aoa.map(c => c[1]))
    ]
}
function drawSmoothPath(ctx, points, rate = 0.2) {
    const add = (p1, p2) => ([p1[0] + p2[0], p1[1] + p2[1]])
    const subtract = (p1, p2) => ([p1[0] - p2[0], p1[1] - p2[1]])
    const curveTo = (p1, p2, p3) => { ctx.bezierCurveTo(p1[0], p1[1], p2[0], p2[1], p3[0], p3[1]) }
    let cp = [] //control point
    for (let i = 0; i < points.length - 2; i++) {
        let p = subtract(points[i + 2], points[i])
        cp.push([Math.round(p[0] * rate), Math.round(p[1] * rate)])
    }
    points.map((p, i, ps) => {
        switch (i) {
            case 0:
                ctx.lineTo(p[0], p[1]); break;
            case 1:
                curveTo(ps[0], subtract(p, cp[i - 1]), p); break;
            case ps.length - 1:
                curveTo(add(ps[i - 1], cp[i - 2]), p, p); break;
            default:
                curveTo(add(ps[i - 1], cp[i - 2]), subtract(p, cp[i - 1]), p); break;
        }
    })
}
const EasingFunctions = {
    // no easing, no acceleration
    linear: t => t,
    // accelerating from zero velocity
    easeInQuad: t => t * t,
    // decelerating to zero velocity
    easeOutQuad: t => t * (2 - t),
    // acceleration until halfway, then deceleration
    easeInOutQuad: t => t < .5 ? 2 * t * t : -1 + (4 - 2 * t) * t,
    // accelerating from zero velocity 
    easeInCubic: t => t * t * t,
    // decelerating to zero velocity 
    easeOutCubic: t => (--t) * t * t + 1,
    // acceleration until halfway, then deceleration 
    easeInOutCubic: t => t < .5 ? 4 * t * t * t : (t - 1) * (2 * t - 2) * (2 * t - 2) + 1,
    // accelerating from zero velocity 
    easeInQuart: t => t * t * t * t,
    // decelerating to zero velocity 
    easeOutQuart: t => 1 - (--t) * t * t * t,
    // acceleration until halfway, then deceleration
    easeInOutQuart: t => t < .5 ? 8 * t * t * t * t : 1 - 8 * (--t) * t * t * t,
    // accelerating from zero velocity
    easeInQuint: t => t * t * t * t * t,
    // decelerating to zero velocity
    easeOutQuint: t => 1 + (--t) * t * t * t * t,
    // acceleration until halfway, then deceleration 
    easeInOutQuint: t => t < .5 ? 16 * t * t * t * t * t : 1 + 16 * (--t) * t * t * t * t
}
function getShiftedCoords(coords, width) {
    var x2, y2, x42, y42
    var shift_coord = []
    for (let i = 0; i < coords.length; i++) {
        if (i == 0) {
            var y1 = coords[i][0]
            var x1 = coords[i][1]
            y2 = coords[i + 1][0]
            x2 = coords[i + 1][1]
            var s124 = width / Math.sqrt(Math.pow(x2 - x1, 2) + Math.pow(y2 - y1, 2))
            x42 = (y2 - y1) * s124
            y42 = -(x2 - x1) * s124
            shift_coord.push([y1 + y42, x1 + x42])
        } else if (i == coords.length - 1) {
            shift_coord.push([y2 + y42, x2 + x42])
        } else {
            var y3 = coords[i + 1][0]
            var x3 = coords[i + 1][1]
            var s235 = width / Math.sqrt(Math.pow(x3 - x2, 2) + Math.pow(y3 - y2, 2))
            var x52 = (y3 - y2) * s235
            var y52 = -(x3 - x2) * s235
            var radian = Math.atan2(x42 * y52 - y42 * x52, x42 * x52 + y42 * y52)
            var scalar = 1 / 2 / Math.pow(Math.cos(radian / 2), 2)
            var sx = (x42 + x52) * scalar
            var sy = (y42 + y52) * scalar
            shift_coord.push([y2 + sy, x2 + sx])
            y2 = coords[i + 1][0]
            x2 = coords[i + 1][1]
            x42 = x52
            y42 = y52
        }
    }
    return shift_coord
}
function perpendicular(p, p1, p2) {
    var x = p1[0],
        y = p1[1],
        dx = p2[0] - p1[0],
        dy = p2[1] - p1[1]
    if (dx !== 0 || dy !== 0) {
        var t = ((p[0] - x) * dx + (p[1] - y) * dy) / (dx * dx + dy * dy) //scale of projection
        if (t > 1) {
            x = p2[0]
            y = p2[1]
        } else if (t > 0) {
            x += dx * t
            y += dy * t
        }
    }
    dx = p[0] - x
    dy = p[1] - y
    return {
        distance: Math.sqrt(dx * dx + dy * dy),
        closest: [x, y]
    }
}
function radialDistFilter(points, epsilon) {
        if (epsilon == 0) return points
        var epsilon_squ = epsilon * epsilon
        return points.filter((point, i, points) => {
            if (i == 0) return true
            var dx = point[0] - points[i - 1][0]
            var dy = point[1] - points[i - 1][1]
            return dx * dx + dy * dy > epsilon_squ
        })
    }
function DouglasPeucker(points, epsilon) {
    // points = this.radialDistFilter(points,epsilon)
    if (epsilon == 0) return points
    var dmax = 0, index = -1, firstPoint = points[0], lastPoint = points[points.length - 1]
    for (var i = 1; i < points.length; i++) {
        var d = this.perpendicular(points[i], firstPoint, lastPoint).distance
        if (d > dmax) {
            dmax = d
            index = i
        }
    }
    if (dmax >= epsilon) {
        return [
            ...this.DouglasPeucker(points.slice(0, index + 1), epsilon).slice(0, -1),
            ...this.DouglasPeucker(points.slice(index), epsilon)
        ]
    } else {
        return [points[0], points[points.length - 1]]
    }
}
```
## Geolocation 網頁定位 API
```graphviz
digraph{
navigator->geolocation;
geolocation->getCurrentPosition,watchPosition,clearWatch;
watchPosition->fn_success,fn_error,options;
fn_success->position;
position->coords;
coords->latitude,longitude,accuracy,altitude,heading,speed,timestamp;
options->enableHighAccuracy,timeout,maximumAge;
}
```