TypeScript + Node.js + D3.js を使って白地図に自転車で訪れた市町村のみを塗りつぶすプログラムを作った

はじめに

自分は自転車でツーリングで行った後は 思い出にふけるために白地図に対して通った市町村などをポチポチ塗りつぶしたりしてニヤニヤしたりしてました。

↓ ポチポチしてたときに使ったサイト

n.freemap.jp

↓ ニヤニヤしてたときのツイート

通過した市町村

千歳市
恵庭市
札幌市
小樽市
余市町
赤井川村
倶知安町
ニセコ町
蘭越町
黒松内町
長万部町
今金町
せたな町
八雲町
森町
七飯町
函館市
北斗市
豊浦町
洞爺湖町
伊達市
室蘭市
登別市
白老町
苫小牧市

赤が宿泊、緑が通過 pic.twitter.com/KXxI2CFHoK

— うっひょい/サイモン・FujiwaLaTeX (@keisuke495500) 2022年5月8日

Strava APIを使ってツーリングのデータを取り出せるのでこれも自動化できそうだなって思ってプログラムを作ってみました。

何を使うか

最初はGoogle Chart APIというデータを図にしたりできるAPIの一つであるGoogle GeoChartを使おうと考えました。

Google GeoChartは日本地図に対して色を塗ったりできますが、塗りつぶしの最小単位が都道府県で、 市町村を塗りつぶすという要件は満たしませんでした。

更に調べていくとD3.jsというJavaScriptのパッケージを使って日本地図を描画している人がちらほらいました。

hiyuzawa.jp

qiita.com

記事も充実してやりやすそうだと思い今回はD3.jsを使うことにしました。

D3.jsとは

D3.jsはデータを視覚化したりすることに使われることが多いですが、D3.js自体に描画機能はなく JSONやCSVなどによって与えたデータから描画位置などを計算したり、その描画位置などのデータを付加してDOM操作ができるライブラリです*1。

描画自体はHTML,CSS,SVGなどのWeb標準の機能に準拠します。

d3js.org

D3.jsには緯度経度から描画用のピクセル座標に変換するメソッドなどもあり、今回はそのあたりのメソッドをメインに使っていきます。

D3.jsで白抜き地図を出力してみる

D3.jsを使って白抜きの地図を出力してみます。そのためには元となる地図データが必要です。

地図データの取得

今回の地図データはROIS-DS人文学オープンデータ共同利用センターが提供している市区町村及び行政区毎に境界線が引いてある北海道の地図データを使用します。

geoshape.ex.nii.ac.jp

↑を見るとどのような地図データなのかのプレビューが見られます。

複数の年月日で地図データがおいてあります。 一番古い1920-01-01と最新の2021-01-01を見比べると市町村の合併などで行政区の境界線が変わっていることがわかります。 今回は最新の2021-01-01を使います。

各年月日でも解像度毎にデータセットが分かれています。あまりにも大きい解像度のデータを使うとこの後の変換処理などが重くなるので、今回は中解像度を使います。

ダウンロードしたファイルの形式がtopojsonというテキストファイルになっています。 この形式のデータを加工してD3.jsでマッピングできる形に持っていきます。

TopoJSONとGeoJsonとは

ダウンロードしたtopojsonという形式のファイルはGISデータを表現できるTopoJSONという規格で保存されたファイルです。 同じくGISデータを表現するためのGeoJSONという規格から軽量化を考えて作られた規格です。

github.com

TopoJSONを読み込む

ダウンロードしたTopoJSONを読み込みます。 TypeScriptで扱うための型定義があるのでインストールします。

www.npmjs.com

TopoJSONはその名の通り、フォーマットはJSONベースなのでJSONファイルとして読み込めます。 ただし、読み込む際の型はTopology型にします。

import { Topology } from "topojson-specification";

async function main() {
  const topo: Topology = JSON.parse(
    fs.readFileSync("./topojson/01_city.i.topojson", "utf-8")
  );
  console.log(topo);
}

読み込んだtopojsonを出力すると、オブジェクトになっていることがわかります。

SVGを出力する

TopoJSONのデータから地図をSVG形式で表示します。

TopoJSONからGeoJSONに変換する

D3.jsで読み込むためにはTopoJSONからGeoJSONに変換する必要があります。 変換にはtopojson-clientライブラリを使うため型定義と共にインストールします。

www.npmjs.com

www.npmjs.com

topojson-clientのfeatureメソッドでTopoJSONの指定したオブジェクトをGeoJSONに変換することができます。 TopoJSONにはGeometryCollectionという要素があります。

qiita.com

GeometryObjectと呼ばれる点、線、面、ポリゴン、GeometryCollection*2などのデータが集まったものがGeometryCollectionです。

GeoJSONにてGeometryCollectionに対応するのがFeatureCollection、GeometryObjectに対応するのがFeature及びFeatureCollectionです。

featureメソッドの第一引数にはtopojsonそのものを、第二引数にはGeometryObjectを指定します。 返り値は指定したGeometryObjectから変換されたFeatureCollectionまたはFeatureです。

今回のtopojsonはobjects.cityにGeometryCollectionが入っています。 試しに出力してみます。

import { Topology } from "topojson-specification";

async function main() {
  const topo: Topology = JSON.parse(
    fs.readFileSync("./topojson/01_city.i.topojson", "utf-8")
  );
  console.log(topo.objects.city);
}

すると、以下のようなオブジェクトになっているのがわかります。

{
  type: 'GeometryCollection',
  geometries: [
    {
      type: 'MultiPolygon',
      arcs: [Array],
      id: '北海道札幌市中央区',
...

featureメソッドを使ってGeometryCollectionをGeoJSONのFeatureCollectionに変換します。

import * as topojson from "topojson-client";
import { Topology } from "topojson-specification";

async function main() {
  const topo: Topology = JSON.parse(
    fs.readFileSync("./topojson/01_city.i.topojson", "utf-8")
  );
  const geo = topojson.feature(topo, topo.objects.city);
  console.log(geo);
}

これで出力すると

{
  type: 'FeatureCollection',
  features: [
    {
      type: 'Feature',
      id: '北海道札幌市中央区',
      properties: [Object],
      geometry: [Object]
    },
    {
...

変換できていることがわかります。

FeatureCollectionはFeatureの集まりです。D3.jsで変換処理する際はFeatureを一つずつ取り出して処理するためFeatureCollectionからFeatureのIterableな値を取り出す必要があります。

FeatureCollectionが持つfeaturesというプロパティにFeatureの配列が入っています。

つまり、geo.featuresでFeatureの配列として取り出せます。

const features = geo.features

しかし、TypeScriptでは、topojson.featureによって変換した際のgeoの型はFeature | FeatureCollectionとなっています*3。 Feature型だった場合はgeo.featuresというプロパティは存在しません。

Feature型かFeatureCollection型かどうかはgeo.typeにstring型で入っています。

よって、featuresは

const features =
  geo.type === "FeatureCollection"
    ? geo.features
    : geo.type === "Feature"
      ? [geo]
      : undefined;
if (!features) {
  throw new Error();
}

typeの場合分けによって取得できます。

GeoJSONを元にSVGを出力する

SVGはXMLによって記述されています。D3.jsでピクセル座標に変換後はDOM操作でSVGを構築する必要がありますが、 今回はNode.jsを使用しているため、Node.jsでDOM操作ができるJSDOMをインストールします。

www.npmjs.com

JSDOMとD3.jsを使ってGeoJSONからSVGに変換します。

DOM操作するためにJSDOMのインスタンスを生成してdocumentオブジェクトを取得します。

const document = new JSDOM().window.document;

次に緯度経度から特定のピクセル座標に変換する関数をD3.jsで作ります。

import * as d3 from "d3";
...
const aProjection = d3
  .geoMercator()
  .center([center.lng, center.lat])
  .translate([width / 2, height / 2])
  .scale(scale);

geoMercator関数は緯度経度をメルカトル図法の地図にマッピングします。 マッピングする際に以下の情報を与えることでその情報を踏まえてマッピングする関数を返します。

• center
  • 中心となるピクセル座標(translate)に対応する緯度経度
• translate
  • 中心となるピクセル座標
• scale
  • 緯度経度からピクセル座標にマッピングする際の倍率

例えば、800x800の画像に5000倍の大きさでSVGを出力したい場合を考えます。

800x800の画像の中心は(400,400)なのでtranslateにはそれを与えます。

5000倍の大きさなのでscaleもそのまま与えます。

centerについて、まずは北海道の地図を出力するので北海道の中心を見つけます。

Google Mapsでだいたい北海道の真ん中あたりに合わせてURLに入っている緯度経度をcenterに与えます。

今回の場合は43.4259796,142.6960534です。 コードにすると以下です。

import * as d3 from "d3";
...
const aProjection = d3
  .geoMercator()
  .center([142.6960534, 43.4259796])
  .translate([400, 400)
  .scale(5000);

centerに与える際は経度、緯度の順番です。 Google Mapsに表示されている順番と逆であることに注意してください。

これで緯度経度からピクセル座標に変換するaProjectionという関数を取得できます。

次にD3.jsのgeoPath関数を使って、GeoJSONのFeatureからSVGのフォーマットに変換する関数を取得します。

geoPath().projection() の引数に先程取得したaProjectionを渡すと、FeatureからSVGフォーマットに変換する関数を取得できます。

import * as d3 from "d3";
...
const aProjection = d3
  .geoMercator()
  .center([142.6960534, 43.4259796])
  .translate([400, 400)
  .scale(5000);

const geoPath = d3.geoPath().projection(aProjection);

geoPathという関数を取得しました。

geoPathをGeoJSONのFeature一つ一つに適用していきSVGフォーマットに変換していきます。 また、同時にDOM操作によってSVGを構築していきます。

  const d3Svg = d3
    .select(document.body)
    .append("svg")
    .attr("xmlns", "http://www.w3.org/2000/svg")
    .attr("width", 800)
    .attr("height", 800);

d3.selectでタグを選択できます。bodyの中にSVGを作りたいので、document.bodyを選択します。

appendでタグを追加できます。svgタグを追加します。

更にattrでsvgタグに属性を追加します。svgタグに必要な属性を追加します。

これで、800x800の空のSVGが出来上がります。

elementオブジェクトとしてd3Svgを取得し、GeoJSONのデータを変換しつつ入れていきます。

  d3Svg
    .selectAll("path")
    .data(features)
    .enter()
    .append("path")
    .attr("d", geoPath)
    .style("stroke", "#000000")
    .style("stroke-width", 3)
    .style("fill", "#ffffff");

dataとenterでpathタグの要素にGeoJSONを入れていきます。

selectAllとappendでpathタグをfeaturesの個数だけ生成します。

attrにgeoPath関数を指定して、featuresの中身一つ一つにgeoPath関数を適用した結果をpathの中に入れていきます。

この辺りの挙動は難しく

wizardace.com

ここを参考にしました。

styleでpathタグにスタイルを適用します。SVGの線の色と太さを指定しています。

これらのDOM操作によって、SVGが構築できました。

実際にHTMLを表示してみます。

console.log(document.body.innerHTML);

すると、

<svg xmlns="http://www.w3.org/2000/svg" width="800" height="800"><path d="M376.59773479337764,408.5635356562202L376.6120848716737,408.56283232223666L376.7894585138747,
...
</path></svg>

このようにSVGが出力されます。 ファイルに出力してブラウザなどで開くと市町村の境界線付きで北海道の地図が表示されます。

指定の市町村によって塗りつぶす色を変える

指定した市町村の区域を別の色で塗りつぶしてみます。

今回ダウンロードしたTopoJSONにはGeometry毎に市町村名がついています。

GeoJSONに変換した後にはFeature毎に対応した市町村名が入っています。

{
  type: 'FeatureCollection',
  features: [
    {
      type: 'Feature',
      id: '北海道札幌市中央区',
      properties: [Object],
      geometry: [Object]
    },
    {
...

idとして入っています。

このidをSVGを構築する際のHTMLのpathタグのclass名として利用します。

  d3Svg
    .selectAll("path")
    .data(features)
    .enter()
    .append("path")
    .attr("d", geoPath)
    .attr("class", (d) => {
      return d.id ? d.id : "unknown";
    })
    .style("stroke", "#000000")
    .style("stroke-width", 1)
    .style("fill", "#ffffff");

pathタグに classという属性を指定します。attrの第二引数に関数を指定できます。

ここに指定した関数の引数であるdはfeaturesから取り出されたデータです。つまりFeatureです。

今回はclass名に市町村名を使うためidを取り出して返す関数をattrの第二引数に指定します。

これで各pathタグに市町村名のclassが付与されました。 selectを使ってCSSセレクタで特定のクラスを指定します。

今回は北海道北見市を青色に塗りつぶしてみます。

d3Svg.select(".北海道北見市").style("fill", "blue");

これで北見市を塗りつぶせます。

出力すると

北見市が青色に塗りつぶされていることがわかります。

d3Svg.select(".北海道北見市").style("fill", "blue");
d3Svg.select(".北海道宗谷郡猿払村").style("fill", "red");
d3Svg.select(".北海道雨竜郡幌加内町").style("fill", "yellow");
d3Svg.select(".北海道釧路市").style("fill", "green");

複数塗りつぶすことも可能です。

飛び地の釧路市もちゃんと塗りつぶせてます。

StravaのGPSデータから訪問した市町村を取り出して塗りつぶす色を変える

次はStravaのGPSデータから訪問した市町村を取り出して塗りつぶす方法を考えます。

Stravaのデータの取得方法やどのようなデータが入っているかは

muscle-keisuke.hatenablog.com

を見てください。

Stravaのデータから得られる位置情報は緯度経度のみです。特定の緯度経度がどの市町村なのかを取得する必要があります。

geoloniaさんが出している緯度経度から市町村を取り出すOSSがあるのでこちらを利用させてもらいます。

github.com

インストールして使ってみます。緯度経度から室蘭工業大学がどこの市町村にあるかを取得してみます。

室蘭工業大学の緯度経度は42.3785905,141.0373005なので

import { openReverseGeocoder } from "@geolonia/open-reverse-geocoder";

async function main() {
  const res = await openReverseGeocoder([141.0373005, 42.3785905]);
  console.log(res);
}

とすれば取得できます。ここでも緯度経度の順番がGoogle Mapsと関数に渡す順番が違うことに注意してください。

出力は

{ code: '01205', prefecture: '北海道', city: '室蘭市' }

です。室蘭工業大学が室蘭市にあることがわかりました。

これでStravaの位置情報を入れることでどこの市町村を通ったかがわかるようになりました。

今回は今年のGWに北海道をツーリングしたときに走ったある1日のデータを使ってみます。

ちなみにこの日は赤井川村からせたな町まで走りました。

あらかじめAPIで取得していたこの日のデータをJSONに保存していたので、これを読み込んで位置情報データを出力してみます。

async function main() {
  const activity = JSON.parse(
    fs.readFileSync("./json/activity.json", "utf-8")
  );
  console.log(activity.map.summary_polyline);

すると、

wqweGksszYtKUdRoa@hGqGhCsUrG...EdKjAxRhPvS_@~LsPbV

文字列が返ってきます。これは位置情報を文字列にして表したPolyline Encodingと呼ばれるものです。 詳しくは過去の記事に載っています。

muscle-keisuke.hatenablog.com

このPolyline Encodingされた位置情報を緯度経度の位置情報に変換する必要があります。 これも変換できるOSSがあるので利用させてもらいます。

github.com

インストールして使います。

import polyline from "@mapbox/polyline";

async function main() {
  const activity = JSON.parse(
    fs.readFileSync("./json/activity.json", "utf-8")
  );

  const positions = polyline.decode(activity.map.summary_polyline);

  console.log(positions);
}

出力すると

[
  [ 43.05196, 140.84422 ], [ 43.04993, 140.84433 ], [ 43.04686, 140.84985 ],
  [ 43.04553, 140.85122 ], [ 43.04484, 140.85484 ], [ 43.04346, 140.85771 ],
  [ 43.04133, 140.85959 ], [ 43.03861, 140.86008 ], [ 43.03689, 140.86253 ],
  [ 43.03551, 140.8635 ],  [ 43.02964, 140.86588 ], [ 43.02341, 140.86626 ],
  [ 43.02118, 140.86784 ], [ 43.01437, 140.86788 ], [ 43.01134, 140.86901 ],
...

緯度経度の配列が取得できました。

一つずつ市町村を取得する関数に通します。

async function main() {
  const activity = JSON.parse(
    fs.readFileSync("./json/activity.json", "utf-8")
  );

  const positions = polyline.decode(activity.map.summary_polyline);

  const regions = await Promise.all(
    positions.map(async (p) => await openReverseGeocoder([p[1], p[0]]))
  );

  console.log(regions);

すると、

[
  { code: '01409', prefecture: '北海道', city: '余市郡赤井川村' },
  { code: '01409', prefecture: '北海道', city: '余市郡赤井川村' },
  { code: '01409', prefecture: '北海道', city: '余市郡赤井川村' },
  { code: '01409', prefecture: '北海道', city: '余市郡赤井川村' },
  { code: '01409', prefecture: '北海道', city: '余市郡赤井川村' },
... 
  { code: '01400', prefecture: '北海道', city: '虻田郡倶知安町' },
  { code: '01400', prefecture: '北海道', city: '虻田郡倶知安町' },
  { code: '01400', prefecture: '北海道', city: '虻田郡倶知安町' },
  { code: '01395', prefecture: '北海道', city: '虻田郡ニセコ町' },
  { code: '01395', prefecture: '北海道', city: '虻田郡ニセコ町' },
  { code: '01395', prefecture: '北海道', city: '虻田郡ニセコ町' },
...
  { code: '01394', prefecture: '北海道', city: '磯谷郡蘭越町' },
  { code: '01394', prefecture: '北海道', city: '磯谷郡蘭越町' },
  { code: '01394', prefecture: '北海道', city: '磯谷郡蘭越町' },
  { code: '01394', prefecture: '北海道', city: '磯谷郡蘭越町' },
  { code: '01394', prefecture: '北海道', city: '磯谷郡蘭越町' },

取得できましたが、一つ一つの緯度経度に対して、市町村を返してしまうので、重複があります。 重複を排除します。

また、codeは必要なくprefectureとcityは一つの文字列に結合してしまっても問題ないのでその処理も合わせて行います。

async function main() {
  const activity = JSON.parse(
    fs.readFileSync("./json/activity.json", "utf-8")
  );

  const positions = polyline.decode(activity.map.summary_polyline);

  const regions = await Promise.all(
    positions.map(async (p) => {
      const res = await openReverseGeocoder([p[1], p[0]]);
      return res.prefecture + res.city;
    })
  );

  const distinctRegions = [...new Set(regions)];

  console.log(distinctRegions);

出力すると

[
  '北海道余市郡赤井川村',
  '北海道虻田郡倶知安町',
  '北海道虻田郡ニセコ町',
  '北海道磯谷郡蘭越町',
  '北海道寿都郡黒松内町',
  '北海道山越郡長万部町',
  '北海道瀬棚郡今金町',
  '北海道久遠郡せたな町'
]

重複を削除して必要な形でデータを取得することができました。

あとは前のセクションで説明した塗りつぶしのコードと組み合わせれば、GPSのデータから訪問した市町村を塗りつぶした地図を出力することができます。

import { Topology } from "topojson-specification";
import * as topojson from "topojson-client";
import * as d3 from "d3";
import { JSDOM } from "jsdom";
import fs from "fs";
import { openReverseGeocoder } from "@geolonia/open-reverse-geocoder";
import polyline from "@mapbox/polyline";

async function main() {
  const document = new JSDOM().window.document;
  const topo: Topology = JSON.parse(
    fs.readFileSync("./topojson/01_city.i.topojson", "utf-8")
  );
  const geo = topojson.feature(topo, topo.objects.city);

  const features =
    geo.type === "FeatureCollection"
      ? geo.features
      : geo.type === "Feature"
        ? [geo]
        : undefined;
  if (!features) {
    throw new Error();
  }

  const aProjection = d3
    .geoMercator()
    .center([142.6960534, 43.4259796])
    .translate([400, 400])
    .scale(5000);
  const geoPath = d3.geoPath().projection(aProjection);
  const d3Svg = d3
    .select(document.body)
    .append("svg")
    .attr("xmlns", "http://www.w3.org/2000/svg")
    .attr("width", 800)
    .attr("height", 800);

  d3Svg
    .selectAll("path")
    .data(features)
    .enter()
    .append("path")
    .attr("d", geoPath)
    .attr("class", (d) => {
      return d.id ? d.id : "unknown";
    })
    .style("stroke", "#000000")
    .style("stroke-width", 1)
    .style("fill", "#ffffff");

  const activity = JSON.parse(
    fs.readFileSync("./json/activity.json", "utf-8")
  );

  const positions = polyline.decode(activity.map.summary_polyline);

  const regions = await Promise.all(
    positions.map(async (p) => {
      const res = await openReverseGeocoder([p[1], p[0]]);
      return res.prefecture + res.city;
    })
  );

  const distinctRegions = [...new Set(regions)];

  for (const region of distinctRegions) {
    d3Svg.select(`.${region}`).style("fill", "#5EAFC6");
  }

  fs.writeFileSync("./hoge.svg", document.body.innerHTML);

出力されたSVGを見てみます。

GPS通りの市町村が塗りつぶされています。

中心の緯度経度をTopoJSONから算出する

緯度経度からピクセル座標に変換する際の中心はGoogle Mapsで大体の中心を調べて手入力してました。

  const aProjection = d3
    .geoMercator()
    .center([142.6960534, 43.4259796])
    .translate([400, 400])
    .scale(5000);

これをTopoJSONから算出します。TopoJSONの中にはbboxというプロパティがあります。今回使用しているTopoJSONの場合は

{"type":"Topology","id":"x0401:01","metadata":{"type":["行政区境界"],"dc:title":"北海道 市区町村:中解像度TopoJSON","dc:source":"N03-21_01_210101.shp","dc:issued":"2021-01-01",..."bbox":[139.33396016902668,41.351645558995415,148.89440319085463,45.55724341395245],
...

このように入っています。 取り出すときは

  const topo: Topology = JSON.parse(
    fs.readFileSync("./topojson/01_city.i.topojson", "utf-8")
  );
  console.log(topo.bbox);

これで取り出せます。取り出した値を出力すると

[
  139.33396016902668,
  41.351645558995415,
  148.89440319085463,
  45.55724341395245
]

緯度経度のペア2つが出力されます。この緯度経度は今回使用している北海道のTopoJSONのバウンディングボックスです。

e-words.jp

北海道の南西端の緯度経度と北東端の緯度経度が入っています。

Google Mapsで出力してみると端と端であることがわかります。

これを利用して中心の緯度経度を求めます。

2点間の緯度経度から中間点を求める方法を解説しているサイトがあります。

tma.main.jp

地球を半径1の球体と考えて緯度経度を極座標として考えます。

詳しい導出方法は元のサイトに書いてあるのでここには書かないですが、手順としては

• 2点の緯度経度をそれぞれ直交座標に変換する
• 変換した直交座標の中間点を求める
• 中間点を緯度経度に変換する

です。

まずは2点の緯度経度を直交座標に変換します。

import { Topology } from "topojson-specification";

async function main() {
  const topo: Topology = JSON.parse(
    fs.readFileSync("./topojson/01_city.i.topojson", "utf-8")
  );

  if (!topo.bbox) {
    throw new Error();
  }

  console.log(latLngToPixels({ lat: topo.bbox[1], lng: topo.bbox[0] }));
  console.log(latLngToPixels({ lat: topo.bbox[3], lng: topo.bbox[2] }));
}

function latLngToPixels(arg: { lat: number; lng: number }) {
  const [radLat, radLng] = [degreeToRadian(arg.lat), degreeToRadian(arg.lng)];
  return {
    x: Math.cos(radLng) * Math.cos(radLat),
    y: Math.cos(radLat) * Math.sin(radLng),
    z: Math.sin(radLat),
  };
}

function degreeToRadian(degree: number) {
  return degree * (Math.PI / 180);
}

Mathの三角関数はラジアンを引数に取るので、変換しています。

結果は

{
  x: -0.5693979186700274,
  y: 0.48917259406333097,
  z: 0.6606785780026414
}
{
  x: -0.5995197259339476,
  y: 0.36173329843931645,
  z: 0.7139503617313299
}

です。

これの中間点を求めます。

const aPixels = latLngToPixels({ lat: topo.bbox[1], lng: topo.bbox[0] });
const bPixels = latLngToPixels({ lat: topo.bbox[3], lng: topo.bbox[2] });
// 正規化するので1/2する必要なし
const midPixels = {
  x: aPixels.x + bPixels.x,
  y: aPixels.y + bPixels.y,
  z: aPixels.z + bPixels.z,
};
const midVectorLength = Math.sqrt(
  midPixels.x ** 2 + midPixels.y ** 2 + midPixels.z ** 2
);
const normalizedMidPixels = {
  x: midPixels.x / midVectorLength,
  y: midPixels.y / midVectorLength,
  z: midPixels.z / midVectorLength,
};
console.log(normalizedMidPixels);

出力すると中間点の直交座標が求まります。

{
  x: -0.5859244396192019,
  y: 0.4265198327143081,
  z: 0.6890380129994881
}

求めた中間点の直交座標を緯度経度に戻します。

...
  const midLat = radianToDegree(Math.asin(normalizedMidPixels.z));
  const midLng = radianToDegree(
    Math.atan2(normalizedMidPixels.y, normalizedMidPixels.x)
  );

  console.log({ midLat, midLng });
}
...
function radianToDegree(radian: number) {
  return (radian * 180) / Math.PI;
}

直交座標は三角関数で求めたので緯度経度に戻すときは逆三角関数を使います。

出力は

{ midLat: 43.55400742272156, midLng: 143.94750122681177 }

これをGoogle Mapsで見てみると、

北方四島を含めると確かに中心に見えます。

求めた中心点をD3.jsに入れてSVG出力してみます。

北方四島含めいバランスよく出力されています。

描画範囲全体がバランスよく出力される倍率を求める

D3.jsに指定するscaleも目視で確認して5000倍という値を入れていますが、これもバウンディングボックスを利用して求めることができます。

800x800のSVGを出力することを考えます。

const width = 800;
const height = 800;

const originProjection = d3.geoMercator().scale(1);

const minPosition = originProjection([bbox[0], bbox[1]]) ?? [0, 0];
const maxPosition = originProjection([bbox[2], bbox[3]]) ?? [0, 0];

const originWidth = Math.abs(maxPosition[0] - minPosition[0]);
const originHeight = Math.abs(maxPosition[1] - minPosition[1]);

const scale = Math.min(width / originWidth, height / originHeight);

まず1倍でメルカトル図法で緯度経度をピクセル座標にマッピングする関数を作ります。

const width = 800;
const height = 800;

const originProjection = d3.geoMercator().scale(1);

そして、バウンディングボックスの2点をピクセル座標に変換します。

const minPosition = originProjection([bbox[0], bbox[1]]) ?? [0, 0];
const maxPosition = originProjection([bbox[2], bbox[3]]) ?? [0, 0];

倍率1倍における端の座標が求まります。端2点のwidthとheightを差から求めます。

const originWidth = Math.abs(maxPosition[0] - minPosition[0]);
const originHeight = Math.abs(maxPosition[1] - minPosition[1]);

必要なSVGの長さは800x800なので、originWidthとoriginHeightの何倍に当たるかを求めます。 これが倍率です。

ただし、倍率そのままだと余白がない状態なので求めたscaleを95%にします。 また、元のバウンディングボックスが縦と横で大きさが異なるため、求まる倍率も異なります。

倍率が大きい方に合わせると当然描画範囲をはみ出るのでMath.minで小さい方を取得します。

const scale = 0.95 * Math.min(width / originWidth, height / originHeight);

北海道は北方四島の分があり、縦よりも横に長くなるため、倍率は横の方が小さくなります。

求めたscaleを使ってSVGを出力してみます。

const aProjection = d3
  .geoMercator()
  .center([midLng, midLat])
  .translate([width / 2, height / 2])
  .scale(scale);

全体がバランス良く出力されています。

SVGからPNGに変換する

Twitterに載せたりするために、SVGからPNGに変換します。 Node.jsで画像処理できるsharpを使ってみます。

www.npmjs.com

型定義も入れます。

www.npmjs.com

出力したSVGファイルに対して、変換を行います。

import sharp from "sharp";

async function main() {
...
  fs.writeFileSync("./hoge.svg", document.body.innerHTML);

  await sharp("./hoge.svg").png().toFile("./hoge.png");
}

800x800のPNG画像が出力されました。

*1:D3.jsはData-Driven-Documentsから名付けられている

*2:GeometryCollectionの入れ子がありえる

*3:featureメソッドにはGeometryObjectが引数として渡るため