🙆♀️
GeoJSONのbboxの座標変換
GeoJSON（EPSG:4326）のbboxを別の座標系に計算し直したいとき、単純に左下と右上の座標をそれぞれ投影してしまいがちですが、これでは正しい座標にならないときがあります。
{
type: "Polygon",
coordinates: [
[
[-180, 80],
[180, 80],
[180, 90],
[-180, 90],
[-180, 80]
]
]
}
例えば極を囲むようなポリゴン[1]のbboxをEPSG:3995に変換したいの場合、4326では[-180, 80, 180, 90]ですが、これを単純に変換してしまうと四角が潰れてしまいます（実際にやってみよう）。
また線分が長いと変換後の線分が曲がり（大圏）、元の四隅からはみ出る部分が発生するかもしれない[2]。
じゃあどうすればいいのかというと、愚直に線分間に点をたくさん挿入して、それらを投影変換後、最大最小を求めてやるくらいしか方法がない。[3]
以下のサンプルコードではGeoJSONの座標は必ずEPSG:4326の緯度経度座標系であると仮定しています。 [4]
function toBbox(
polygon,
code
) {
const xs = [];
const ys = [];
if (code && code !== "EPSG:4326") {
for (let i = 0; i < polygon.coordinates[0].length - 1; i++) {
const interpolated = linearInterpolationPoints(
polygon.coordinates[0][i].slice(0, 2),
polygon.coordinates[0][i + 1].slice(0, 2),
{
partition: 10, // エイヤの数字
}
);
for (let j = 0; j < interpolated.length - 1; j++) {
const transformed = proj4("EPSG:4326", code, [
interpolated[j][0],
interpolated[j][1],
]) as [number, number];
xs.push(transformed[0]);
ys.push(transformed[1]);
}
}
} else {
polygon.coordinates[0].forEach((lonlats) => {
xs.push(lonlats[0]);
ys.push(lonlats[1]);
});
}
const minLon = Math.min(...xs);
const minLat = Math.min(...ys);
const maxLon = Math.max(...xs);
const maxLat = Math.max(...ys);
return [minLon, minLat, maxLon, maxLat];
}
function calcLinearInterpolationPoint (
i,
p1,
p2,
partition
) {
const div = Math.floor(partition);
if (_eq(p1[0], p2[0])) {
return [p1[0], p1[1] + (i * (p2[1] - p1[1])) / (div + 1)];
} else if (_eq(p1[1], p2[1])) {
return [p1[0] + (i * (p2[0] - p1[0])) / (div + 1), p1[1]];
}
const x = p1[0] + (i * (p2[0] - p1[0])) / (div + 1);
return [x, linearInterpolationY(p1, p2, x)];
};
function linearInterpolationPoints(
p1,
p2,
userOptions={}
): Points {
const options = Object.assign(
{
partition: 0,
},
userOptions
);
if (options.partition <= 0) return [p1, p2];
const ret: Points = [];
for (let i = 0; i < options.partition + 1; i++) {
ret.push(calcLinearInterpolationPoint(i, p1, p2, options.partition));
}
ret.push(p2);
return ret;
}
function linearInterpolationY(
p1,
p2,
x
) {
if (_eq(p1[0], p2[0])) return p1[1];
return p1[1] + ((x - p1[0]) * (p2[1] - p1[1])) / (p2[0] - p1[0]);
}
function _eq(a, b) {
return Math.abs(a - b) < 1e-8;
}
やってることは基本的にはこの記事でやってることの逆の操作。
toBbox(
{
type: "Polygon",
coordinates: [
[
[-180, 80],
[180, 80],
[180, 90],
[-180, 90],
[-180, 80],
],
],
},
"EPSG:3995"
)
> [-1078093.1792349985, -1045060.0346785864, 1078093.1792349985, 1089179.4556261837]
Discussion