地球上の夜の領域をGeoJSON化する
その時刻の夜(または昼)の領域を世界地図上に表示したくなる瞬間が人生にはある。
この昼夜の境界線を英語ではTerminatorという[1]。
境界線の座標は幾何学的に求めることができる。
計算方法は日の出・日の入りの計算―天体の出没時刻の求め方に詳しい[2]。
ウェブページではこちらがよくまとまっている。
出没方程式もしくは薄明方程式
太陽が地平線の下に隠れてもしばらく明るい薄明時間というが、
地平線に太陽が隠れた瞬間を求めたければ
なお、
よって
解析解
となり解析解が求まる。
数値解
時角
太陽の赤経・赤緯
赤経は時角を求めるために必要。
順序としては、
- 計算したい日時からユリウス通日を求める
- ユリウス通日から太陽の黄経を求める
- ユリウス通日から黄道傾斜角を求める
- 太陽の黄経と黄道傾斜角から赤経・赤緯を求める
各式は「日の出・日の入りの計算―天体の出没時刻の求め方」やwikipediaを参照のこと。
プログラム
import { Feature, Polygon } from "geojson";
const getJulianDate = (date: Date | number): number => {
// https://www5d.biglobe.ne.jp/~noocyte/Programming/GregorianAndJulianCalendars.html
const d = date instanceof Date ? date.getTime() : date;
return d / 86400000 + 2440587.5;
}
const getGreenwichMeanSiderealTime = (jd: number): number => {
// https://aa.usno.navy.mil/faq/GAST
return (18.697374558 + 24.06570982441908 * (jd - 2451545)) % 24;
}
const getSunEclipticLongitude = (jd: number): number => {
// https://en.wikipedia.org/wiki/Position_of_the_Sun
const n = jd - 2451545;
const L = (280.46 + 0.9856474 * n) % 360;
const g = ((357.528 + 0.9856003 * n) % 360) * Math.PI / 180;
// radians
return (L + 1.915 * Math.sin(g) + 0.02 * Math.sin(2 * g)) * Math.PI / 180;
}
const getEclipticObliquity = (jd: number): number => {
// https://ja.wikipedia.org/wiki/%E9%BB%84%E9%81%93%E5%82%BE%E6%96%9C%E8%A7%92
const t = (jd - 2451545) / 36525;
// radians
return ((84381.406 - 46.836769 * t - 0.00059 * t ** 2 + 0.001813 * t ** 3) / 3600) * Math.PI / 180;
}
const getSunEquatorialPosition = (longitude: number, obliquity: number): number[] => {
// 日の出・日の入りの計算 天体の出没時刻の求め方 長沢工著 ISBN4-8052-0634-9
let alpha =
Math.atan(Math.tan(longitude) * Math.cos(obliquity));
const delta =
Math.asin(Math.sin(longitude) * Math.sin(obliquity));
if (Math.sin(longitude) > 0) {
if (Math.sin(alpha) < 0) alpha += Math.PI
} else {
if (0 < Math.sin(alpha)) alpha += Math.PI
}
// radians
return [alpha, delta];
}
const getLatitude = (longitude: number, alpha: number, delta: number, gmst: number): number => {
// 日の出・日の入りの計算 天体の出没時刻の求め方 長沢工著 ISBN4-8052-0634-9
// 太陽高度kを0に固定
const hourAngle = (gmst * 15 * Math.PI / 180 + longitude) - alpha;
// radians
return Math.atan(-Math.cos(hourAngle) / Math.tan(delta));
}
const calculate = (date: Date, division = 360): Feature<Polygon> => {
const jd = getJulianDate(date);
const gmst = getGreenwichMeanSiderealTime(jd);
const sunEclipticLongitude = getSunEclipticLongitude(jd);
const eclipticObliquity = getEclipticObliquity(jd);
const [alpha, delta] = getSunEquatorialPosition(sunEclipticLongitude, eclipticObliquity);
const longlats: number[][] = [];
const diffDeg = 360 / division;
for (let longitude = -180; longitude < 180; longitude += diffDeg) {
const latitude = getLatitude(longitude * Math.PI / 180, alpha, delta, gmst);
longlats.push([longitude, latitude * 180 / Math.PI]);
}
const latitude = getLatitude(Math.PI, alpha, delta, gmst);
longlats.push([180, latitude * 180 / Math.PI]);
if (delta > 0) {
longlats.push([180, -90]);
longlats.push([-180, -90]);
longlats.push(longlats[0]);
longlats.reverse();
} else {
longlats.push([180, 90]);
longlats.push([-180, 90]);
longlats.push(longlats[0]);
}
return {
type: "Feature",
properties: {
"datetime": date.toISOString()
},
geometry: {
type: "Polygon",
coordinates: [longlats]
}
}
}
GeoJSONのPolygonはcounterclockwiseである必要があるので南半球側に陰があるときはlonglats.reverse();
している
2023-10-01T07:00:00Z
ごろのポリゴンをgeojson.ioで表示してみた様子
-
検索時は某映画と区別がつくよう「terminator day night line」などと言葉を足すとよい。 ↩︎
-
長沢工先生の本はすべて良著で、これだけの内容を日本語で学べることに感謝しきりである。 ↩︎
-
薄明時間の分類は世界的にはこの三種類が一般的で、
は一等星が見えだす限界の明るさ、k=-6 は航海中に水平線が見える限界の明るさ、k=-12 は六等星が見えだす限界の明るさということらしい。 ↩︎k=-18 -
本来は線上の地点の標高や視差も方程式に含まれるが、ここではどちらも0としている。 ↩︎
-
「日の出・日の入りの計算―天体の出没時刻の求め方」のP58では緯度を2deg刻みで計算すれば十分とされている。 ↩︎
Discussion