-
Notifications
You must be signed in to change notification settings - Fork 9
/
02_merge.php
122 lines (118 loc) · 4.25 KB
/
02_merge.php
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
<?php
$jsonPath = __DIR__ . '/geojson';
if(!file_exists($jsonPath)) {
mkdir($jsonPath, 0777);
}
foreach(glob(__DIR__ . '/raw/*.json') AS $jsonFile) {
echo "processing {$jsonFile}\n";
$json = json_decode(file_get_contents($jsonFile), true);
// if(!isset($json['features'])) {
// echo "{$jsonFile}\n";
// print_r($json); exit();
// }
// continue;
foreach($json['features'] AS $f) {
$targetPath = $jsonPath . '/' . substr($f['attributes']['段號'], 0, 6);
if(!file_exists($targetPath)) {
mkdir($targetPath, 0777);
}
$targetFile = $targetPath . '/' . $f['attributes']['段號'] . '.json';
$feature = array(
'type' => 'Feature',
'properties' => array(
'段號' => $f['attributes']['段號'],
'地籍址' => $f['attributes']['地籍址'],
'PLNInt' => $f['attributes']['PLNInt'],
'農盤分類' => $f['attributes']['農盤分類'],
'使用分區' => $f['attributes']['使用分區'],
'使用類別' => $f['attributes']['使用類別'],
'Shape_Length' => $f['attributes']['Shape_Length'],
'Shape_Area' => $f['attributes']['Shape_Area'],
),
'geometry' => array(),
);
// if(!isset($f['geometry']['rings'])) {
// print_r($json); exit();
// }
if(count($f['geometry']['rings']) === 1) {
$feature['geometry'] = array(
'type' => 'Polygon',
'coordinates' => array(),
);
foreach($f['geometry']['rings'] AS $k1 => $ring) {
foreach($ring AS $k2 => $point) {
$point = twd97_to_latlng($point[0], $point[1]);
$ring[$k2][0] = floatval($point['lng']);
$ring[$k2][1] = floatval($point['lat']);
}
$feature['geometry']['coordinates'][] = $ring;
}
} else {
$feature['geometry'] = array(
'type' => 'MultiPolygon',
'coordinates' => array(),
);
foreach($f['geometry']['rings'] AS $k1 => $ring) {
foreach($ring AS $k2 => $point) {
$point = twd97_to_latlng($point[0], $point[1]);
$ring[$k2][0] = floatval($point['lng']);
$ring[$k2][1] = floatval($point['lat']);
}
$feature['geometry']['coordinates'][] = array($ring);
}
}
file_put_contents($targetFile, json_encode($feature));
}
}
foreach(glob(__DIR__ . '/geojson/*') AS $tmpPath) {
$fc = array(
'type' => 'FeatureCollection',
'features' => array(),
);
foreach(glob($tmpPath . '/*.json') AS $jsonFile) {
$fc['features'][] = json_decode(file_get_contents($jsonFile));
unlink($jsonFile);
}
file_put_contents($tmpPath . '.json', json_encode($fc));
rmdir($tmpPath);
}
function twd97_to_latlng($x, $y) {
$a = 6378137.0;
$b = 6356752.314245;
$lng0 = 121 * M_PI / 180;
$k0 = 0.9999;
$dx = 250000;
$dy = 0;
$e = pow((1 - pow($b, 2) / pow($a, 2)), 0.5);
$x -= $dx;
$y -= $dy;
$M = $y / $k0;
$mu = $M / ($a * (1.0 - pow($e, 2) / 4.0 - 3 * pow($e, 4) / 64.0 - 5 * pow($e, 6) / 256.0));
$e1 = (1.0 - pow((1.0 - pow($e, 2)), 0.5)) / (1.0 + pow((1.0 - pow($e, 2)), 0.5));
$J1 = (3 * $e1 / 2 - 27 * pow($e1, 3) / 32.0);
$J2 = (21 * pow($e1, 2) / 16 - 55 * pow($e1, 4) / 32.0);
$J3 = (151 * pow($e1, 3) / 96.0);
$J4 = (1097 * pow($e1, 4) / 512.0);
$fp = $mu + $J1 * sin(2 * $mu) + $J2 * sin(4 * $mu) + $J3 * sin(6 * $mu) + $J4 * sin(8 * $mu);
$e2 = pow(($e * $a / $b), 2);
$C1 = pow($e2 * cos($fp), 2);
$T1 = pow(tan($fp), 2);
$R1 = $a * (1 - pow($e, 2)) / pow((1 - pow($e, 2) * pow(sin($fp), 2)), (3.0 / 2.0));
$N1 = $a / pow((1 - pow($e, 2) * pow(sin($fp), 2)), 0.5);
$D = $x / ($N1 * $k0);
$Q1 = $N1 * tan($fp) / $R1;
$Q2 = (pow($D, 2) / 2.0);
$Q3 = (5 + 3 * $T1 + 10 * $C1 - 4 * pow($C1, 2) - 9 * $e2) * pow($D, 4) / 24.0;
$Q4 = (61 + 90 * $T1 + 298 * $C1 + 45 * pow($T1, 2) - 3 * pow($C1, 2) - 252 * $e2) * pow($D, 6) / 720.0;
$lat = $fp - $Q1 * ($Q2 - $Q3 + $Q4);
$Q5 = $D;
$Q6 = (1 + 2 * $T1 + $C1) * pow($D, 3) / 6;
$Q7 = (5 - 2 * $C1 + 28 * $T1 - 3 * pow($C1, 2) + 8 * $e2 + 24 * pow($T1, 2)) * pow($D, 5) / 120.0;
$lng = $lng0 + ($Q5 - $Q6 + $Q7) / cos($fp);
$lat = ($lat * 180) / M_PI;
$lng = ($lng * 180) / M_PI;
return array(
'lat' => round($lat, 7),
'lng' => round($lng, 7)
);
}