1.はじめに

私が出没する「QGIS初心者質問グループ」に国土基本図図郭を作成するサイトが動作しなくなっているとの投稿がありました。確かにEdge、Chrome、Firefoxのいずれでも、投稿のとおりプルダウンができず、動作しません。またググってみたところ、容易に図郭を入手できるような手段もみつかりませんでした。
そこでQGISを使っての手作業で国土基本図図郭を作成し、投稿に添付しましたが、作成手順について説明したものがあった方が良いと思い、これを機会に記事にまとめることにしました。
なお国土基本図図郭に関する説明は、空間情報クラブさんの情報がとてもわかりやすかったです。私の説明でわかりにくいところがあれば、そちらをご覧頂ければと思います。

2.国土基本図郭について

国土基本図図郭は、平面直角座標を元に定義されています。平面直角座標は原点を持ちますので、原点からの距離で長方形の枡に区切り、それぞれの枡をルールに従ってネーミングします。
国土基本図図郭には地図情報レベルが5段階（50000、5000、2500、1000、500）あるようですが、ここでは地図情報レベル50000、5000、2500の3種類について説明します。
なお図郭のネーミングルールでは冒頭には系番号が付くことになっていますが、この記事では省略しています。

2-1.地図情報レベル50000

国土基本図図郭のエリアは、平面直角座標で用いられる19の各原点から南北へそれぞれ300km、東西にそれぞれ160kmです。
これを南北方向30km、東西方向40kmの枠で区切ったものが「地図情報レベル50000」です。
北から南へ「A～T」、西から東へ「A～H」が割り当てられ、図郭コードは「南北アルファベット＋東西アルファベット」となります。

2-2.地図情報レベル5000

地図情報レベル5000は、地図情報レベル50000の図郭を東西南北方向に10等分にしたものです。北から南へ「0～9」、西から東へ「0～9」が割り当てられ、図郭コードは「地図情報50000図郭コード」＋「南北番号＋東西番号」となります。

2-3.地図情報レベル2500

地図情報レベル2500は、地図情報レベル5000の図角を東西南北方向に2等分したものです。北西「1」、北東「2」、南西「3」、南東「4」が割り振られ、図郭コードは「地図情報5000図郭コード」の後ろにこの番号を付加したものになります.

3.図郭ポリゴンの生成

ここからが作業手順になります。QGISを立ち上げましょう。
QGISが立ち上がったら、最初にプロジェクトのプロパティでCRSを図角を生成したい地域の平面直角座標に設定します。
次の例では、JGD2011の平面直角座標10系に設定しています。

JGD2000とJGD2011の平面直角座標のEPSGは次表のとおりです。プロジェクトのCRSを設定する際などの参考としてください。

系番号	JGD2000 (EPSG)	JGD2011 (EPSG)	適用区域
1	2443	6669	長崎県鹿児島県のうち北方北緯32度南方北緯27度西方東経128度18分東方東経130度を境界線とする区域内（奄美群島は東経130度13分までを含む。)にあるすべての島、小島、環礁及び岩礁
2	2444	6670	福岡県　佐賀県　熊本県　大分県　宮崎県　鹿児島県（I系に規定する区域を除く。)
3	2445	6671	山口県　島根県　広島県
4	2446	6672	香川県　愛媛県　徳島県　高知県
5	2447	6673	兵庫県　鳥取県　岡山県
6	2448	6674	京都府　大阪府　福井県　滋賀県　三重県　奈良県和歌山県
7	2449	6675	石川県　富山県　岐阜県　愛知県
8	2450	6676	新潟県　長野県　山梨県　静岡県
9	2451	6677	東京都（14系、18系及び19系に規定する区域を除く。)　福島県　栃木県　茨城県　埼玉県千葉県　群馬県　神奈川県
10	2452	6678	青森県　秋田県　山形県　岩手県　宮城県
11	2453	6679	小樽市　函館市　伊達市　北斗市　北海道後志総合振興局の所管区域　北海道胆振総合振興局の所管区域のうち豊浦町、壮瞥町及び洞爺湖町　北海道渡島総合振興局の所管区域　北海道檜山振興局の所管区域
12	2454	6680	北海道（11系及び13系に規定する区域を除く。）
13	2455	6681	北見市　帯広市　釧路市　網走市　根室市　北海道オホーツク総合振興局の所管区域のうち美幌町、津別町、斜里町、清里町、小清水町、訓子府町、置戸町、佐呂間町及び大空町　北海道十勝総合振興局の所管区域　北海道釧路総合振興局の所管区域　北海道根室振興局の所管区域
14	2456	6682	東京都のうち北緯28度から南であり、かつ東経140度30分から東であり東経143度から西である区域
15	2457	6683	沖縄県のうち東経126度から東であり、かつ東経130度から西である区域
16	2458	6684	沖縄県のうち東経126度から西である区域
17	2459	6685	沖縄県のうち東経130度から東である区域
18	2460	6686	東京都のうち北緯28度から南であり、かつ東経140度30分から西である区域
19	2461	6687	東京都のうち北緯28度から南であり、かつ東経143度から東である区域

3-1.地図情報レベル50000の図郭作成

メニューの「ベクタ－調査ツール－グリッドを作成」でポリゴンを生成します。
次図でグリッドの範囲は「-160000,160000,-300000,300000 [EPSG:6678]」としていますが、生成する地域に合わせてEPSGのコード番号は、適宜修正してください。
また、「水平方向の間隔」は「40000」、「垂直方向の間隔」は「30000」としています。

生成されたデータを確認してみます。
属性テーブルをみると、ID番号と上下左右の座標値が入っています。

フィールド計算機を使って座標値からコードを生成しますが、考え方を順に追って説明します。

まず行ですが、上から下へ数が増えていくように番号を振りたいです。
ここでは座標値を30000で除して得た値を使えそうですが、正の値となるところでは上から下へ数が減っていきます。そこで「top」というフィールドから300000を減じた値を30000で除したものの絶対値を使うことにします。これで、きれいに上から下へ向かって0～19という数値が振られます。

次に列も左から右へ数が増えて行くように振りたいです。
さきほどの考えを応用して、「left」というフィールドへ160000を加算した値を40000で除した値を使います。これで、左から右へ向かって0～7という数値が振られます。

続いて数値を英文字にする方法ですが、文字列関数の「char」を使います。これはユニコード番号に対応した文字を変えす関数で、「A」のユニコード番号「65」です。
よって割り振った数値に65を加算したすると、目的の英文字に対応するユニコード番号が得られます。

以上のことを踏まえると、フィールド計算機でコードを求める式は次のようになります。

char( abs( "top" - 300000) / 30000 + 65 ) || char(( "left" + 160000) / 40000 + 65 )

3-2.地図情報レベル5000の図郭作成

地図情報レベル50000の時と同様の方法でポリゴンを作成しますが、「水平方向の間隔」は「4000」、「垂直方向の間隔」は「3000」とします。なお、次の属性の空間結合を行うという工程でファイルを出力しますので、「グリッド」という項目の設定はデフォルトの「一時レイヤを作成」のままにします。

ポリゴンを生成したら、先ほど作成した地図情報レベル50000のコードを属性テーブルに取り込みます。
メニューの「ベクタ－データ管理ツール－属性の空間結合」を立ち上げ、次図のような設定で実行します。

地図情報レベル50000の後に続く番号ですが、縦方向を3000、横方向を4000で割って得た値だと、縦が0～199、横方向が0～79になります。
そこで、さらに10で割った余りの値（剰余）を使うように式を組み立てます。
※「code50000」は空間結合で属性テーブルに取り込んだ地図情報レベル50000のコードが入っているフィールド

"code50000" || to_string(( abs(( "top" - 300000 )/ 3000 ))% 10) || to_string((("left" + 160000)/ 4000) % 10)

3-3.地図情報レベル2500の作成

ポリゴンを「水平方向の間隔」を「2000」、「垂直方向の間隔」を「1500」で作成して、属性の空間結合で地図レベル5000からコードを取得するところまでは同様の手順となります。
ただし地図情報レベル2500のコードは、行番号と列番号を振って組み合わせるのではなく、一連の番号（1～4）を振る必要があるので工夫が必要となります。

まず行番号と列番号を0と1の繰り返しにすることは、先ほどまでと同様の手順でできると思います。
そこでちょっと工夫です。行番号が0と2の繰り返し、列番号が0と1の繰り返しだったら、行番号と列番号を足すと0～3の数値が得られます。これに1を足すと1～4の連番となり、目的の数値が得られます。
以上のことを式に組み立てます。
※「code5000」は空間結合で属性テーブルに取り込んだ地図情報レベル5000のコードが入っているフィールド

"code5000" || to_string( ((abs(( "top" - 300000 )/ 1500 ))% 2)* 2 +(("left" + 160000)/ 2000) % 2 + 1)

4.さいごに

フィールド計算機にコピペして一気に計算できるように、コードを得るための式は1本にまとめることにしました。これで作業工程数は減りますが、内容は少しわかりにくいかもしれません。
このような際、私は以下のように式を分解して出力を確認することを良くやります。

内容を確認する際には、それぞれ計算していると思われる部分をコピペしてフィールド計算機を実行し、出力の様子をラベルで表示して確認すると良いと思います。
例えば私が「地図情報レベル50000」の式を確認する場合ですと、次のようなります。

1.「abs( "top" - 300000) / 30000」で行番号の計算結果について確認
2.「1」の式からabsを取って、絶対値でなければ計算結果がどのようになるかを確認
3.「1」の式で「-300000」が無ければ計算結果がどのようになるかを確認
4.「( "left" + 160000) / 40000」で列番号の計算結果について確認
5.「char("「1」の計算結果のフィールド" + 65)」で行に対応するアルファベットへの変換の確認
6.「char("「4」の計算結果のフィールド" + 65)」で列に対応するアルファベットへの変換の確認

実際の作業では「1」「4」「5」「6」と作業を行って、「5」と「6」を「||」で文字列を繋いで結果を得ても何ら差し支えないと思います。この方が記憶しやすいですし、中途の計算過程を確認できる堅実な方法だと思います。