前回の記事（国土数値情報の令和2年の医療圏データの文字化けが直っていた）では国土数値情報の医療圏データで四苦八苦した様子を掲載した．今回は二次医療圏テーブルの文字化けの解決を試みる．
データの中身は？
相手を知るには SQL Server に用意されているメソッドを使ってプロパティを探る．今回は STGeometryType() メソッドと NumRings() メソッドを使ってみる．前者はジオメトリの種類を返し，後者はポリゴンのリングの総数を返すが，マルチポリゴンだと NULL を返す仕様である．
下記クエリを各テーブルに対して実行する．
SELECT geog.STGeometryType() AS GeometryType , geog.NumRings() AS NumRings FROM [MedicalAreaDB].[dbo].[T_MedicalArea1] SELECT geog.STGeometryType() AS GeometryType , geog.NumRings() AS NumRings FROM [MedicalAreaDB].[dbo].[T_MedicalArea2] SELECT geog.STGeometryType() AS GeometryType , geog.NumRings() AS NumRings FROM [MedicalAreaDB].[dbo].[T_MedicalArea3]
結果は以下の通りである．
(118119 行処理されました) (116365 行処理されました) (116037 行処理されました)
一次医療圏テーブルにマルチポリゴンが一つ使用されているが，残りはすべてポリゴンである．NumRings() メソッドの結果が興味深いため掲載しておく．
一次医療圏
| NumRings()メソッドの返り値 | 個数 |
|---|---|
| 1 | 118018 |
| 2 | 69 |
| 3 | 18 |
| 4 | 5 |
| 5 | 4 |
| 6 | 2 |
| 9 | 1 |
| 11 | 1 |
| NULL | 1 |
二次医療圏
| NumRings()メソッドの返り値 | 個数 |
|---|---|
| 1 | 116336 |
| 2 | 21 |
| 3 | 5 |
| 4 | 1 |
| 5 | 1 |
| 6 | 1 |
三次医療圏
| NumRings()メソッドの返り値 | 個数 |
|---|---|
| 1 | 116025 |
| 2 | 9 |
| 3 | 1 |
| 4 | 1 |
| 6 | 1 |
QGISで融合（ディゾルブ）を実行する
二次医療圏は335件，三次医療圏は52件に集約するのが本来のデータベースの形である．その結果空間データはマルチポリゴンになるであろうことは予想できる．問題は属性テーブルの文字化けである．地物を表現する空間データと属性テーブルを分けてアップロードし，SQL Server 内で結合してはどうだろうか．二次医療圏の属性は下記クエリで抽出できる．
WITH CTE
AS
(
SELECT DISTINCT
[A38b_001]
,[A38b_002]
,[A38b_003]
,[A38b_004]
,[A38b_005]
,[A38b_006]
,[A38b_007]
,[A38b_008]
,[A38b_009]
,[A38b_010]
,[A38b_011]
FROM [MedicalAreaDB].[dbo].[T_MedicalArea2]
)
SELECT * FROM CTE
(335 行処理されました)
前回の記事（国土数値情報の令和2年の医療圏データの文字化けが直っていた）に従って二次医療圏テーブル T_MedicalArea2 のポリゴンを融合しよう．二次医療圏コード (A38b_003) を集約キーとして融合(dissolve)する．
二次医療圏を表現する335件のマルチポリゴンができる．
「エクスポート」「新規ファイルに地物を保存」
下図を参考に融合した地物をエクスポートする．
- 形式を ESRI Shapefile とする．
- ファイル名をフルパスで指定する．ここでは T_MedicalArea2_Dissolved としている．
- 座標参照系を EPSG:4326 に変更する．
- ジオメトリをポリゴンとし，マルチタイプにするをチェックする．
下図を参考に属性テーブルをエクスポートする．
- 形式をカンマで区切られた値[CSV]とする．
- ファイル名をフルパスで指定する．ここでは T_MedicalArea2_TEXT としている．
- 座標参照系を EPSG:4326 とする．
- 文字コードを Shift_JIS とする．
- ジオメトリを「ジオメトリなし」とする．
- レイヤオプションで LINEFORMAT を CRLF とする．
- SEPARATPR を TAB とする．
- STRING_QUOTED を IF_NEEDED にする．
ウィザードを使って属性テーブルをSQL Serverにインポート
- T_MedicalArea2_TEXT をメモ帳で開き，ダブルクォーテーションマークを一括置換で削除する
- T_MedicalArea2_TEXT の拡張子をtxtに変更する
- SQL Server のウィザードを使って T_MedicalArea2_TEXT を SQL Server にインポートする
Shape2SQLでSQL Serverにアップロード
下図を参考に Shape2SQL で SQL Server にアップロードする．
- Geography を選択する
- Set SRID をチェックして 4326 を指定する
ジオメトリの種類と数を調べる
STGeometryType() メソッドはジオメトリの種類を返す．下記クエリを実行してみる．
SELECT [geog].STGeometryType() AS GeometryType , LEN(geog.STAsText()) AS LENGTH , geog.NumRings() AS NumRings FROM [MedicalAreaDB].[dbo].[T_MedicalArea2_Dissolved]
結果をExcelにコピペして GeometryType の個数をピボットテーブルに表示する．
| GeometryType | 個数 |
|---|---|
| GeometryCollection | 2 |
| MultiPolygon | 204 |
| Polygon | 129 |
一次医療圏のポリゴンを二次医療圏コードで集約し融合した場合，GeometryCollection の数が 3 に増え，MultiPolygon の数が 203 に減ったが Polygon の数は 129 で変わらなかった．
NumRings() メソッドは Polygon に適用された場合に限って数値を返すが，MultiPolygon に適用されても NULL を返す仕様である．下表は NumRings の集計結果である．
| NumRings | 個数 |
|---|---|
| 1 | 123 |
| 2 | 6 |
| NULL | 206 |
STAsText() メソッドはジオメトリの WKT 表現を返す．その文字列長を返す LENGTH の結果は割愛する．
UPDATE SET FROM JOIN ステートメントでテーブルを更新する
この構文は初めて使う．手始めに次のクエリを実行して結果を確認する．
SELECT M2.geog , T.* FROM dbo.T_MedicalArea2_Dissolved AS M2 INNER JOIN dbo.T_MedicalArea2_Text AS T ON M2.A38b_003 = T.A38b_003
結果は以下の通り．
(335 行処理されました)
更新クエリ
もともと同じレイヤからエクスポートしたデータであり，二次医療圏コードで結合して元のテーブルに戻したいのである．次のクエリを実行する．
USE MedicalAreaDB; GO UPDATE M2 SET M2.A38b_001 = T.A38b_001 , M2.A38b_002 = T.A38b_002 , M2.A38b_004 = T.A38b_004 , M2.A38b_005 = T.A38b_005 , M2.A38b_006 = T.A38b_006 , M2.A38b_007 = T.A38b_007 , M2.A38b_008 = T.A38b_008 , M2.A38b_009 = T.A38b_009 , M2.A38b_010 = T.A38b_010 , M2.A38b_011 = T.A38b_011 FROM dbo.T_MedicalArea2_Dissolved AS M2 INNER JOIN dbo.T_MedicalArea2_Text AS T ON M2.A38b_003 = T.A38b_003
(335 行処理されました)
下記クエリを実行して結果を確認する．
SELECT TOP 10 M2.A38b_001 , M2.A38b_002 , M2.A38b_003 , M2.A38b_004 FROM dbo.T_MedicalArea2_Dissolved AS M2
A38b_001 A38b_002 A38b_003 A38b_004 35207,35210,35215 下松市,光市,周南市 3503 周南 35203,35206 山口市,防府市 3504 山口・防府 35208,35321 岩国市,和木町 3501 岩国 22211,22213,22216,22223,22224,22461 磐田市,掛川市,袋井市,御前崎市,菊川市,森町 2207 中東遠 35212,35305,35341,35343,35344 柳井市,周防大島町,上関町,田布施町,平生町 3502 柳井 22131,22132,22133,22134,22135,22136,22137,22221 中区,東区,西区,南区,北区,浜北区,天竜区,湖西市 2208 西部 35211 長門市 3507 長門 35204,35502 萩市,阿武町 3508 萩 35202,35213,35216 宇部市,美祢市,山陽小野田市 3505 宇部・小野田 35201 下関市 3506 下関
どうやらうまく行ったようである．
結果
ジオメトリを空間結果タブで確認する
例として北海道の二次医療圏を抽出するクエリを実行する．
SELECT * FROM dbo.T_MedicalArea2_Dissolved AS M2 WHERE LEFT(M2.A38b_003, 2) = '01'
結果は下図の通り．ボンヌ図法で表現してあるが，空間結果タブには「1つ以上の空間オブジェクトが大きすぎて表示できませんでした．クエリを修正してください」と警告が表示され，根室医療圏が表示できていない．QGIS なら遅いながらも全て表示可能なのだと思われる．
根室医療圏のデータを深掘りしてみよう．
SELECT geog.STGeometryType() AS GeometryType , geog.STNumPoints() AS NumPoints , geog.STLength() AS Length , geog.STArea() AS Area , geog.STNumGeometries() AS NumGeometries FROM dbo.T_MedicalArea2_Dissolved AS M WHERE M.A38b_003 = '0121'
結果は以下の通りである．
GeometryType NumPoints Length Area NumGeometries MultiPolygon 106850 1008832.61170782 3592178524.40058 1469
統計を取ってみる
下記クエリを実行してジオメトリの種類，Pointの数，周囲長，面積，ジオメトリ数の統計を取ってみよう．
SELECT M.A38b_003 , M.A38b_004 , geog.STGeometryType() AS GeometryType , geog.STNumPoints() AS NumPoints , geog.STLength() AS Length , geog.STArea() AS Area , geog.STNumGeometries() AS NumGeometries FROM dbo.T_MedicalArea2_Dissolved AS M ORDER BY M.A38b_003
結果の例としてジオメトリ数とポイント数の分布を両対数の散布図に示す．先の投稿で問題を提起したが，二次医療圏データは一次医療圏のジオメトリを流用したわけではなく，日本という国土の複雑さがもたらした結果であることの傍証と言えるだろう．
二次医療圏ごとの病床数の合計を調べる
先の記事（国土数値情報の医療機関のPointデータをダウンロードしQGIS経由でSQL Serverにアップロードする）で医療機関の空間データをアップロードした．下記クエリは所属する二次医療圏ごとの各医療機関の病床数の合計を抽出するものである．
SELECT SUM(CAST(M.[P04_008] AS bigint)) AS '病床数合計' , A.A38b_003 AS '二次医療圏コード' , A.A38b_004 AS '二次医療圏名' , A.A38b_007 AS '医療計画人口' , A.A38b_008 AS '総人口' , A.A38b_009 AS '15歳未満人口' , A.A38b_010 AS '15才以上65才未満人口' , A.A38b_011 AS '65歳以上人口' FROM [MedicalAreaDB].[dbo].[T_Medical_Institution] AS M INNER JOIN dbo.T_MedicalArea2_Dissolved AS A ON M.Location.STWithin(A.geog) = 1 GROUP BY A.A38b_003 , A.A38b_004 , A.A38b_007 , A.A38b_008 , A.A38b_009 , A.A38b_010 , A.A38b_011
クエリの実行には恐ろしく時間がかかる（筆者の環境で 43 分）．これは空間演算の STWithin() メソッドが 18 万行×335 行の合計 60,740,860 件を処理しなければならないためである．空間演算は計算コストが高く，しないで済むならそれに越したことはない．国土交通省におかれては，医療機関の属性にせめて二次医療圏コードを付与しておいてもらいたい．
二次医療圏ごとの総人口と病床数の合計とを散布図に示す．総じて人口と病床数には線形の相関関係が認められる．決定係数は 0.8287 と高い．ちなみに15歳未満人口の決定係数は 0.8205, 15歳以上65歳未満人口の決定係数は 0.8057, 65歳以上人口の決定係数は 0.8656 であった．病床数は 65 歳以上人口を基準に決定されている可能性がある．
近似曲線より上のポイント，つまり人口に比して病床数が多い医療圏では供給過多となっている可能性があり，個々の医療機関の経営が苦しいと予想される．逆に近似曲線より下のポイント，つまり人口に比して病床数が少ない医療圏では常に病床が逼迫しており，供給不足となっている可能性がある．いずれにしても医療機関の経営者にとっては頭の痛い問題である．
医療提供体制には病床数の他，医師数もその要因として挙げられる．国土数値情報にはそのデータはなく，別のデータベースを当たる必要がある（例：地域医療情報システムなど）．将来の医療需要などの予測は国際医療福祉大学の高橋らが公開している．
まとめ
二次医療圏テーブルの文字化けを解決するため，ジオメトリと属性に分けてそれぞれを SQL Server にアップロードし，SQL Server 内で属性を更新した．ジオメトリは Shape2SQL を使ってアップロードし，属性はウィザードを使ってインポートした．
先の投稿で提起した疑問，二次医療圏のデータは一次医療圏のデータを流用したのではないかというものだが，そうではなく，日本の国土の複雑さが二次医療圏のデータ件数の大きさの原因であると推測された．
二次医療圏コードで 335 件に集約した結果，空間データはマルチポリゴンまたはジオメトリコレクションの形を取り，一行あたりのポリゴン数の最大値は 7596 個（石巻・登米・気仙沼），ポイント数の最大値は 531130 個（石巻・登米・気仙沼）になった．
“国土数値情報の二次医療圏テーブルの文字化けを解決する” への2件の返信