# GIS - 地圖開發
> GIS 地理資訊系統(Geographic Information System)
## 概念
向量(vector) | 網格(raster)
-|-
縮放不失真、儲存使用量小 | 渲染效率高，容易切割資源檔
### 座標參考系統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;
}
```
{"metaMigratedAt":"2023-06-14T22:59:07.444Z","metaMigratedFrom":"YAML","title":"GIS - 地圖開發","breaks":true,"contributors":"[{\"id\":\"794553c4-57e8-424f-b9da-1f0d40a61e1b\",\"add\":47822,\"del\":27767}]"}